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ABSTRACT 

We report on a method, PUSH, for artificially triggering core-collapse supernova explosions of massive stars 
in spherical symmetry. We explore basic explosion properties and calibrate PUSH to reproduce SN 1987A 
observables. Our simulations are based on the GR hydrodynamics code AGILE combined with the neutrino 
transport scheme IDSA for electron neutrinos and ASL for the heavy flavor neutrinos. To trigger explosions in 
the otherwise non-exploding simulations, the PUSH method increases the energy deposition in the gain region 
proportionally to the heavy flavor neutrino fluxes. We explore the progenitor range 18-21 Mq. Our studies 
reveal a distinction between high compactness (HC) (compactness parameter ^1.75 > 0.45) and low compactness 
(LC) (^1.75 < 0.45) progenitor models, where EC models tend to explode earlier, with a lower explosion energy, 
and with a lower remnant mass. HC models are needed to obtain explosion energies around 1 Bethe, as observed 
for SN 1987A. However, all the models with sufficiently high explosion energy overproduce ^®Ni and fallback 
is needed to reproduce the observed nucleosynthesis yields. ^^“^^Ni yields depend sensitively on the electron 
fraction and on the location of the mass cut with respect to the shell structure of the progenitor. We identify 
a progenitor and a suitable set of parameters that fit the explosion properties of SN 1987A assuming 0.1 Mq 
of fallback. We predict a neutron star with a gravitational mass of 1.50 Mq. We And correlations between 
explosion properties and the compactness of the progenitor model in the explored mass range. However, a 
more complete analysis will require exploring of a larger set of progenitors. 

Subject headings: hydrodynamics, nucleosynthesis, stars: neutron, supernovae; general, supernovae: individ¬ 
ual (SN 1987A) 


1. INTRODUCTION 

Core-collapse supernovae (CCSN) occur at the end of the 
life of massive stars (M > 8 - 10 Mq). In these violent 
events, the core of the star gravitationally collapses and trig¬ 
gers a shock wave, leading to the supernova explosion. De¬ 
spite many decades of theoretical and numerical modeling, 
the detailed explosion mechanism is not yet fully under¬ 
stood. Simulations in spherical symmetry including detailed 
neutrino transport and general relativity fail to explode self- 
consi stently, except for th e lowest-mass core-coll apse progen¬ 
itors (iFischer et al.ll2010t iHiidepohl et al.l 120101) . There are 
many ongoing efforts using multi-dimensional fluid dynam¬ 
ics, magnetic fields, and rotation to address various remaining 
open questio n s in core-co l lapse s upernova theory (see, e.g., 
iJankal (l2012h : iJanka et ^ (120121) : iBurrowsI (l2013h ). Among 
those are also technical issues, for example the consequences 
of neutrino transport approximations, the convergence of sim¬ 
ulation results, or the dependence of the simulation outcome 
on the dimensionality of the model. Awareness of this depen¬ 
dence is especially important because not all investigations 
can be performed in a computationally very expensive three- 
dimensional model. While sophisticated multi-dimensional 
models are needed for an accurate investigation of the ex¬ 
plosion mechanism, they are currently too expensive for sys¬ 
tematic studies that have to be based on a large number of 
progenitor models. But such a large number of simulations 
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is required to address the following fundamental questions; 
What are the conditions for explosive nucleosynthesis as a 
function of progenitor properties? What is the connection be¬ 
tween the progenitor model and the compact remnant? How 
do these aspects relate to the explosion dynamics and ener¬ 
getics? The lack of readily calculable supernova simulations 
with self-consistent explosions is a problem for many related 
fields, in particular for predicting nucleosynthetic yields of 
supernovae. As we will continue to argue below, spherically 
symmetric models of the explosion of massive stars are still 
a pragmatic method to study large numbers of stellar progen¬ 
itors, from the onset of the explosion up to several seconds 
after core bounce. 

In the past, supernova nucleosynthesis predictions re¬ 
lied on artiflcially triggered explosions, either using a pis¬ 
ton (e.g. IWooslev & Weaveijll995l iLimongi & Chieffil 120061 
Chieffi & Limongil l2013h or a thermal energy b omb (e.g., 
Thielemann et al.l 119961 : lUmeda & Nomotol 120081) . For the 
piston model, the motion of a mass point is specified along 
a ballistic trajectory. For the thermal energy bomb, explo¬ 
sions are triggered by adding thermal energy to a mass zone. 
In both cases, additional energy is added to the system to trig¬ 
ger an explosion. In addition, the mass cut (bifurcation be¬ 
tween the proto-neutron star (PNS) and the ejecta) and the 
explosion energy are free parameters which haye to be con¬ 
strained from the mass of the ^®Ni ejecta. While these ap¬ 
proaches are appropriate for the outer layers, where the nu¬ 
cleosynthesis mainly depends on the strength of the shock 
waye, they are clearly incorrect for the innermost layers. 
There, the conditions and the nucleosynthesis are directly re¬ 
lated to the physics of collapse and bounce, and to the de¬ 
tails of the explosion mechanism. Besides the piston and 
thermal bomb methods, another widely used way to artifl- 
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cially trigger explosions is the so-called “neutrino light-bulb”. 
In this method, the PNS is excised and replaced with an 
inner boundary condition which contains an analytical pre¬ 
scription for the neutrino luminosities. The neutrino trans¬ 
port is replaced by neutrino absorption and emission terms 
in optically thin conditions. Suitable choices of the neu¬ 
trino luminosi ties and energies can trigger neutrino-driven ex¬ 
plosions ('e.g..lBur rows & Goshv 1 19 93t I Yamasaki & Yamadal 


l2005tllwakami et al.l2008ll2009l: Yamamoto et al.l2013l) . The 

light-bulb method has also been used to investigate mod¬ 
els with respect to their dimensionality. The transition from 
spherical symmetry (ID) to axisymmetry (2D) delivers the 
new degree of freedom to bring cold accreting matter down to 
the neutrinospheres while matter in other directions can dwell 
longer in the gain re gion and efficiently be heated by neutrinos 


longer in me gam re gion ana emciently be neatea by neutrino; 
(iHerantet alJll994l). The standing accretion shock instabil 
ity (SASI, e.g., Blondin et a l.ll2003l iBlondin (fe'M 


2006: IScheck et ah l2008 r llwakami et al.l 120091 : Fernandez 


ezzacagga 


201C : lGuilet & Foglizzol2012^ strongly contributes to this ef 


feet in 2D light-bulb models and leads to strong polar os- 
cillations of expansion dur ing the unfolding of the explosion 
(iMurphv & Burrowsl2008h . It was first expected that the trend 
toward a smaller critical luminosity for successful explosions 
will con tinue as one goes from 2D to three-dimen sional (3D) 
models (iNordhaus et alJl2010t iHandv et al.l 1201 4), bu t other 
studie s pointed toward the contrary (iHanke et al.l2012HCouchl 
120131) . One has to keep in mind, that a light bulb approach 
might not include the full coupling between the accretion rate 
and the neutrino luminosity. However, recent models that de¬ 
rive the neutrino luminosity from a consistent evolution of the 
neutron star support the res ult that 2D models show faster ex¬ 
plosions than 3D models jMuller et al.||201 ^ iBruenn et alJ 
1201 3t iDolence et al.ll2013l : iTakiwaki et al.ll2014f) . Most im- 


portant for this work is a finding that is consistent with all 
above investigations; In 3D there is no preferred axis. The 3D 
degrees of freedom lead to a more efficient cascade of fluid 
instabilities to smaller scales. In spite of vivid fluid instabil¬ 
ities, the 3D models show in their overall evolution a more 
pronounced sphericity than the 2D models. Hence their av¬ 
erage conditions resemble more closely the shock expansion 
that would be obtained by an exploding ID model. 

In a ID model with detailed Boltzmann neutrino transport 
two other methods to trigger explosions using neutri nos have 
been used (iFrohlich et al.Tl200« iFischer et al.ll20ldl) . These 
“absorption methods” aim at increasing the neutrino energy 
deposition in the heating region by mimicking the expected 
net effects of multi-dimensional simulations. In one case, the 
neutral-current scattering opacities on free nucleons are artifi¬ 
cially decreased to values between 0.1 and 0.7 times the origi¬ 
nal values. This leads to increased diffusive neutrino fluxes in 
regions of very high density. The net results are a faster delep- 
tonization of the PNS and higher neutrino luminosities in the 
heating region. In the other case, explosions are enforced by 
multiplying the reaction rates for neutrino absorption on free 
nucleons by a constant factor. To preserve detailed balance, 
the emission rates also have to be multiplied by the same fac¬ 
tor. This reduces the timescale for neutrino heating and again 
results in a more efficient energy deposition in the heating re¬ 
gion. However, the energy associated with these explosions 
were alway s weak. _ 

Recently, lUgliano et akl (12012h have presented a more so¬ 
phisticated light-bulb method to explode spherically symmet¬ 
ric models using neutrino energy deposition in post-shock lay¬ 
ers. They use an approximate, grey neutrino transport and re¬ 


place the innermost 1.1 M© of the PNS by an inner boundary. 
The evolution of the neutrino boundary luminosity is based 
on an analytic cooling model of the PNS, which depends on 
a set of free parameters. These parameters are set by fitting 
observational propertie s of SN 1987A f or progenitor masses 
around 20 M© (see also lErtl et alJ (1201 Sh i. 

Artificial supernova explosions have been obtained by other 
authors using a grey leakage scheme that includes neutrino 
heating via a parametriz ed charged-current absorption scheme 
( O|Connor_&OdK01T ) in spherically symmetric simulations 
( O’Connor &Ot)ir2011 1. 

In this paper, we report on a new approach, PUSH, for ar¬ 
tificially triggering explosions of massive stars in spherical 
symmetry. In PUSH, we deposit a fraction of the luminos¬ 
ity of the heavy flavor neutrinos emitted by the PNS in the 
gain region to increase the neutrino heating efficiency. We 
ensure an accurate treatment of the electron fraction of the 
ejecta through a spectral neutrino transport scheme for and 
Vg and a detailed evolution of the PNS. We calibrate our new 
method by comparing the explosion energies and nucleosyn¬ 
thesis yields of different progenitor stars with observations of 
SN 1987A. This method provides a framework to study many 
important aspects of core-collapse supernovae for large sets of 
progenitors; explosive supernova nucleosynthesis, neutron- 
star remnant masses, explosion energies, and other aspects 
where full multi-dimensional simulations are still too expen¬ 
sive and traditional piston or thermal bomb models do not cap¬ 
ture all the relevant physics. With PUSH we can investigate 
general tendencies and perform systematic parameter vari¬ 
ations, providing complementary information to “ab-initio” 
multi-dimensional simulations. 

The article is organized as follows; Section |2] describes 
our simulation framework, the stellar progenitor models, the 
new method PUSH, and our post-processing analysis. In Sec- 
tion[3 we present a detailed exploration of the PUSH method 
and the results of fitting it to observables of SN 1987A. We 
also analyze aspects of the supernova dynamics and progeni¬ 
tor dependency. In Section [H we discuss further implications 
of our results and also compare with other works from the lit¬ 
erature. A summary is given and conclusions are drawn in 
Section|5] 


2 . METHOD AND INPUT 
2.1. Hydrodynamics and neutrino transport 

We make use of the general r elativistic hydrodynamics code 
AGILE in spherical symmetry (iLiebendorfer et a02001h . Eor 
the stellar co ll apse, we apply the deleptonization scheme of 
iLiebendorfeii (l2005h . Eor the neutrino transport, we em¬ 
ploy the Isotropic Diffusion Source Approximation (IDSA) 
for the electron neutrino s Vg and electron anti-neutrinos Vg 
(ILiebendorfer et ai1l2009h . and an Advanced Spectral Leak¬ 
age scheme (ASL) for the heavy -lepton flavor neutrinos Vx - 
v^,v^,Vt,Vt (iPerego et al.ll2014l) . We discretize the neutrino 
energy using 20 geometrically increasing energy bins, in the 
range 3MeV < Ey < 300 MeV. The neutrino reactions 
included in the IDSA and ASL scheme are summarized in 
Table [U They represent the minimal set of the most rele¬ 
vant weak processes in the post-bounce phase, particularly 
up to the onset of an explosion. Note that electron cap¬ 
tures on heavy nuclei and neutrino scattering on electrons, 
whic h are relevant in the collapse phase (see, for exam¬ 
ple, iMezzacappa & Bruennlll9933lbl) . are not included ex¬ 
plicitly in the form of reaction rates, but as part of the pa- 
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rameterized deleptonization scheme. In the ASL scheme, we 
omit nucleon-nucleon bremsstrahlung, N + N N + N + 
Vr + Vr, where N denotes any nucleon (se e, for example, 
iHannestad & RafF^ll998l iBartl et aTI 120141) . We have ob¬ 
served that its inclusion would overestimate fj. and t neutrino 
luminosities during the PNS cooling phase, due to the miss¬ 
ing neutrino thermalization provided by inelastic scattering on 
electrons and positrons at the PNS surface. However, we have 
also tested that the omission of this process does not signifi¬ 
cantly change the /i and t neutrino luminosities predicted by 
the ASL scheme before the explosion sets in. Thus, neglect¬ 
ing N-N bremsstrahlung is only relevant for the cooling phase, 
where it improves the overall behavior when compared to sim- 
ulations obtained w ith detailed Boltzmann neutrino transport 
(iFischer et al.ll2010[ e.g.). 

The _ equation _ of _ state (EOS) of 

iHempel & Schaffner-Bielichl (120101) (HS) that we are 
using includes various light nuclei, such as alphas, deuterons 
or tritons (details below). However, the inclusion of all 
neutrino reactions for this detailed nuclear composition 
would be beyond the standard approach implemented in 
current supernova simulations, where only scattering on 
alpha particles is typically included. To not completely 
neglect the contributions of the other light nuclei, we have 
added their mass fractions to the unbound nucleons. This is 
motivated by their very weak binding energies and, therefore, 
by the idea that they behave similarly as the unbound nucleon 
component. 

In all our models we use 180 radial zones which include 
the progenitor star up to the helium shell. This corresponds 
to a radius of R a; (1.3 - 1.5) x 10'° cm. With this setup, we 
model the collapse, bounce, and onset of the explosion. The 
grid of AGILE is adaptive, with more resolution where the 
gradients of the thermodynamic variables are steeper. Thus, in 
the post-bounce and explosion phases, the surface of the PNS 
and the shock are the better resolved regions. The simulations 
are run for a total time of 5 s, corresponding to > 4.6 s after 
core bounce. At this time, the shock has not yet reached the 
external edge of our computational domain. 


Table 1 

Relevant neutrino reactions. 


Reactions 

Treatment 

Reference 

e~ + p <r^ n + Ve 

IDSA 

a 

e* +11 p + Vc 

IDSA 

a 

N + V ++ N + V 

IDSA & ASL 

a 

(A,Z) + v« (A,Z) + v 

IDSA & ASL 

a 

+e* ++ v„,T + Vk.t 

ASL 

a, b 


Note. — Nucle ons are denote d by N. The nucleon charged current 
rates are based on IBruemil il985 |). but effects of mean-field interactions 

IReddv et al.ll99^ : IRoberts et al.l2()12l : IMartmez-Pinedo et al.l2Q12t[HemDell 

are talten into account. 

References. — fai lBrueni] 119851) : (bl lMezzacappa & Bmennl <1993al) . 


2.2. Equation of state and nuclear reactions 

Eor the high-density plasma in nuclear statistical equi¬ 
librium (NSE) the tabulated microphysical EOS HS(DD2) 
is used. This supernova EOS is b ased on the model 
of IHempel & Schaffner-Bielichl (120101) . It use s the DD2 
param etrization for the nucleo n interactions (iTvpel et al.l 
1201 Oh . the nuclear masses fr om lAud i et al.l (1 20031) . and the 
Finite Range Droplet Model (iMoller et al.lll995al) . 8140 nu¬ 
clei are included in total, up to Z = 136 and to the neutron 


drip li ne. The HS(DD2) was first introduced in lEischer et al.l 
(120141) . where its characteristic properties were discussed 
and general EOS effects in core-collapse sup ernova simu¬ 
lations were investigated. lEischer et al.l (1201 4h showed that 
the HS(DD2) EOS gives a better agreement with constraints 
from nuclear experiments and astrophysical observations than 
the commonly us ed EOSs of iLattimer & Swestyl (119911) and 
IShen et aTI (119981) . Eurthermore, additional degrees of free¬ 
dom, such as yarious light nuclei and a statistical ensemble 
of heayy nuclei, are taken into account. The nucleon mean- 
field potentials, which are used in the charged-c urrent rates, 
haye been calculated consistently (lHempelll2014l) . The max¬ 
imum ma ss of a cold neutro n star for the HS(DD2) EOS is 
2.42 Mm (lEisc her et al.l l2014l) . which is well aboye the limits 
from lDemorest et al.l (l2010l) and I Antoni adis et akl (120131) . 

The EOS employed in our simulations includes an exten¬ 
sion to non-NSE conditions. In the non-NSE regime the nu¬ 
clear composition is described by 25 representatiye nuclei 
from neutrons and protons to iron-group nuclei. The cho¬ 
sen nuclei are the alpha-nuclei "^He, '^C, *®0, ^°Ne, ^"^Mg, 
^'^Si, ^^S, ^®Ar, '‘“Ca, '‘^Ti, '^'^Cr, ^^Pe, ^^Ni, complemented 

by '"^N and the following asymmetric isotopes: ^He, ^®S, 
50Ti, 54pg^ 56pg^ 58pg^ 60pg^ 62pg^ 62j^j 

clei it is possible to achieye a mapping of the abundances 
from the progenitor calculations onto our simulations which 
is consistent with the proyided electron fraction, i.e., man- 
taini ng charge neutrality . All the nuclear masses M, are taken 
from lAudi et aTI (120031) . To adyect the nuclear composition 
inside the adaptiye grid, we implem ent the Consisten t Multi - 
fluid Adyection (CMA) method by iPlewa & Miillerl (119991) . 
Eor giyen abundances, the non-NSE EOS is calculated based 
on the same underlying description used in the NSE regime 
(IHempel & Schaffner-Bielichl |2010^ . but with the following 
modifications: excited states of nuclei are neglected, excluded 
yolume effects are not taken into account, and the nucle¬ 
ons are treated as non-interacting Maxwell-Boltzmann gases. 
Such a consistent description of the non-NSE and NSE phases 
preyents spurious effects at the transitions between the two 
regimes. 

Outside of NSE, an approximate a-network is used to fol¬ 
low the changes in composition. Explosiye Helium-, Carbon-, 
Neon-, and Oxygen-burning are currently implemented in the 
simulation. Note that the thermal energy generation by the nu¬ 
clear reactions is fully incorporated yia the detailed non-NSE 
treatment. We do not haye to calculate explicitly any energy 
liberation, but just the changes in the abundances. Within our 
relatiyistic treatment of the EOS (applied both in the non-NSE 
and NSE regime) energy conseryation means that the specific 
internal energy e\ai is not changed by nuclear reactions. This 
is due to fact that e-mt includes the specific rest mass energy 
Cmass, where Cmass IS giycn by the sum oyer the masses of all 
nuclei weighted with their yield Y, = XilAi, 

emass = ^ LMi . (1) 

i 

Howeyer, if we define the specific thermal energy Cth as 

^th ~ ^int “ ^mass7 (2) 

the nuclear reactions will decrease the rest mass energy (i.e., 
increase the binding) and consequently increase the thermal 
energy. This treatment of the non-NSE EOS is consistent with 
the conyention used in all high-density NSE EOSs. 

Due to limitations of our approximate or-network, and be- 





































































4 


Perego et al. 


cause we do not include any quasi-statistical equilibrium de¬ 
scription, we apply some parameterized burning for temper¬ 
atures between 0.3 and 0.4 MeV. A temperature dependent 
burning timescale is introduced, which gradually transforms 
the initial non-NSE composition towards NSE. Eor temper¬ 
atures of 0.4 MeV and above, the non-NSE phase always 
reaches a composition close to NSE and its thermodynamic 
properties become very similar to the NSE phase of the 
HS(DD2) EOS. Even though the two phases are based on the 
same input physics, small, but unavoidable differences can re¬ 
main, due to the limited set of nuclei considered in non-NSE. 
To assure a smooth transition for all conditions, we have in¬ 
troduced a transition region as an additional means of ther¬ 
modynamic stability. We have chosen a parameterization in 
terms of temperature and implement a linear transition in the 
temperature interval from 0.40 MeV to 0.44 MeV. We check 
that the basic thermodynamic stability condition dsjdT > 0 
is always fulfilled. 


2.3. Initial models 

Eor this study, we use solar-metallicity, non-rotating 
stellar models from the stellar evolution code KEPLER 
(IWooslev et al.l |2002). Our set includes 16 pre-explosion 
models with zero-age main sequence (ZAMS) mass between 
18.0 Mq and 21.0 Mq in increments of 0.2 Mq. These models 
have been selected to have ZAMS mass around 20 Mq, sim¬ 
ilar to the progenitor of SN 1987A (e.g.. lPodsiadlowski et al.l 
|2003). We label the models by their ZAMS mass. In Pigure[T] 
the density profiles of the progenitor models are shown. Eor 
each o f them the compactness parameter is defined fol¬ 
lowing lO’ Connor & OttI (1201 ll) by the ratio of a given mass 
M and the radius R{M) which encloses this mass: 


= 


MIMq 

/?(M)/ 1000 km' 


(3) 


Typically, either 75 or ^ 2.5 are used. The compactness can 
be compu ted at the onset of colla pse or at bounce, as sug¬ 
gested by lO’ Connor & Otl (1201 Ih . Eor our progenitors, the 
difference in the compactness parameter between these two 
moments is not significant for our discussions. Thus, for sim¬ 
plicity, in the following we will use f 1 75 computed at the onset 
of the collapse. The progenitor models considered here fall 
into two distinct families of compactness: low compactness 
(^ 1.75 < 0.45; EC models) and high compactness (^ 1.75 > 0.45; 
HC models), see Table |2] Eigure|2] shows the compactness as 
function of ZAMS mass for the progenitors of this study. The 
non-monotonous behavior is a result of the evolution before 
collapse. The mass range between 19 and 21 Mq is particu¬ 
larly prone to variations of the compactness. Eor a detailed 
discussion of th e behavior of the compactnes s as function of 
ZAMS mass see lSukhbold & WooslevI (120141) . 


2.4. The PUSH method 
2.4.1. Rationale 

The goal of PUSH is to provide a computationally efficient 
framework to explode massive stars in spherical symmetry to 
study multiple aspects of core-collapse supernovae. The us¬ 
age of a spectral transport scheme to compute the Ve and 
luminosities provides a more accurate evolution of Yg of the 
innermost ejecta, which is a crucial aspect for nucleosynthe¬ 
sis. The neutrino luminosities include the accretion contribu¬ 
tion, as well as the luminosity coming from PNS mantle and 



Figure 1. Density profiles as function of ZAMS mass for the progenitor 
models included in this study (18.0 Mq to 21.0 M0). HC models are shown 
in red, LC models are shown in blue. The vertical line in the inset is located 
at 1.75 M 0 and indicates that mass at which the compactness parameter ^ 1,75 
is determined (see Equation[3. 
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Figure!. Compactness ^ 1.75 as function of ZAMS mass for our pre¬ 
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core. The accretion luminosity depends not only on the accre¬ 
tion rate but also on the evolution of the mass and radius of 
the PNS, which is treated accurately and self-consistently in 
our models. 

In order to trigger explosions in the otherwise non¬ 
exploding spherically symmetric simulations, we rely on the 
delayed neutrino-driven mechan ism, which was first pro¬ 
posed by iBethe & WilsonI (Il985h . Despite the lack of con¬ 
sensus and convergence of numerical results between dif¬ 
ferent groups, recent multi-dimensional simulations of CC- 
SNe have shown that convection, turbulence and SASI in the 
shocked layers increase the efficiency at which Ve and Vg are 
absorbed inside the gain region, co mpared with spherically 
symmetric models (see, for example. iJanka & Muelleiill996l; 

Nordhat]^_et alJ^10yHanke_e^lj2^2L l2013HDolence et alJ 

2013HCouchll2013HMelson et al.ll2015h . This effect, together 

with the simultaneous increase in time that a fluid particle 
spends inside the g ain region (e.g.. lMurphv & Burrowsll2008t 
iHandv et al.l 120141) . provides more favorable conditions for 
the development of an explosion. Moreover, according to 
multi-dimensional explosion models, the shock revival is fol- 
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Table 2 

Progenitor properties 


Mzams 


.75 

^2.5 


^pr02 

Mpe 

Mco 

Mne 

Mgnv 

(Mo) 

at collapse 

at bounce 

at collapse 

at bounce 

(Mo) 

(Mo) 

(Mo) 

(Mo) 

(Mo) 

18.2 

0.37 

0.380 

0.173 

0.173 

14.58 

1.399 

4.174 

5.395 

9.186 

18.6 

0.365 

0.375 

0.170 

0.170 

14.85 

1.407 

4.317 

5.540 

9.313 

18.8 

0.357 

0.366 

0.166 

0.166 

15.05 

1.399 

4.390 

5.613 

9.435 

19.6 

0.282 

0.288 

0.118 

0.117 

13.37 

1.461 

4.959 

6.243 

7.125 

19.8 

0.334 

0.341 

0.135 

0.135 

14.54 

1.438 

4.867 

6.112 

8.428 

20.0 

0.283 

0.287 

0.125 

0.125 

14.73 

1.456 

4.960 

6.215 

8.517 

20.2 

0.238 

0.241 

0.104 

0.104 

14.47 

1.458 

5.069 

6.342 

8.125 

18.0 

0.463 

0.485 

0.199 

0.199 

14.50 

1.384 

4.104 

5.314 

9.187 

18.4 

0.634 

0.741 

0.185 

0.185 

14.82 

1.490 

4.238 

5.459 

9.366 

19.0 

0.607 

0.715 

0.191 

0.191 

15.03 

1.580 

4.461 

5.693 

9.341 

19.2 

0.633 

0.737 

0.191 

0.192 

15.08 

1.481 

4.545 

5.760 

9.325 

19.4 

0.501 

0.535 

0.185 

0.185 

15.22 

1.367 

4.626 

5.860 

9.365 

20.4 

0.532 

0.594 

0.192 

0.192 

14.81 

1.500 

5.106 

6.376 

8.433 

20.6 

0.742 

0.95 

0.278 

0.279 

14.03 

1.540 

5.260 

6.579 

7.450 

20.8 

0.726 

0.904 

0.271 

0.272 

14.34 

1.528 

5.296 

6.609 

7.735 

21.0 

0.654 

0.764 

0.211 

0.212 

13.00 

1.454 

5.571 

6.969 

6.026 


Note. — ZAMS mass, compactness ^ 1.75 and ^ 2.5 at the onset of collapse and at bounce, total progenitor mass at collapse (Mprog), mass of the iron core (Mpe), 
carbon-oxygen core (Mco). and helium core (Mhc). and mass of the hydrogen-rich envelope (Menv) at collapse, for all the progenitor models included in this 
study. The top part of the table includes the low-compactness progenitors (LC; ^ 1.75 < 0.4 at collapse), the bottom part includes the high-compactness progenitors 
(HC; ^ 1,75 > 0.45 at collapse). 


lowed by a phase where continued accretion a nd shock ex- 


pansion coexist over a time scale of > 1 s (e.e..lScheck et al. 

2006 

;IMarek & Jankal2009HBruenn et al.ll2014UMelson et al. 

2015 

). During this phase, matter accreted through low- 


entropy downflows onto the PNS continues to power an ac¬ 
cretion luminosity. The re-ejection of a fraction of this matter 
by neutrino heating accelerates the shock and increases the 
explosion energy. The length of this phase, the exact amount 
of injected energy, and its deposition rate are still uncertain. 

Inspired by the increase of the net neutrino heating that a 
fluid element experiences due to the above mentioned multi¬ 
dimensional effects, PUSH provides a more efficient neutrino 
energy deposition inside the gain region in spherically sym¬ 
metric models. However, unlike other methods that use elec¬ 
tron flavor neutrinos to trigger artificial ID explosions (see 
Section [TJ, in PUSH we deposit a fraction of the luminos¬ 
ity of the heavy flavor neutrinos (V;f’s) behind the shock to 
ultimately provide successful explosion conditions. This ad¬ 
ditional energy deposition is calibrated by comparing the ex¬ 
plosion energies and nucleosynthesis yields obtained from our 
progenitor sample with observations of SN 1987A. This en¬ 
sures that our artificially increased heating efficiency has an 
empirical foundation. Thus, we can make predictions in the 
sense of an effective model. 

Despite the fact that V;c’s contribute only marginally to the 
energy deposition inside t he gain region in self- consistent 
models (see, for example, iBethe & WilsonI 119851) and that 
they only show a weak dependence on t he temporal varia¬ 
tion o f the accretion rate (see, for example. iLiebendorfer et all 
12004) . their usage presents a number of advantages for our 
purposes. They represent one of the largest energy reser¬ 
voirs available, but they do not directly change the electron 
fraction (unlike electron flavor neutrinos). This allows us 
to trigger an explosion in ID simulations without modify¬ 
ing Ve and Ve luminosities nor changing charged current reac¬ 
tions. The Vx luminosities are calculated consistently within 
our model. They include dynamical feedback from the accre¬ 
tion history, progenitor properties of each individual model. 


and the cooling of the f orming compact object. As shown by 
lO’ Connor & Otd (12013h in broad progenitor studies, during 
the accretion phase that precedes the shock revival, the prop¬ 
erties of the Vx spectral fluxes correlate significantly with the 
properties of v^’s and Vg’s. Unlike the electron (anti-)neutrino 
luminosities, that in spherically symmetric models decrease 
suddenly once the shock has been revived, Vx luminosities are 
only marginally affected by the development of an explosion. 
This allows PUSH to continue injecting energy inside the ex¬ 
panding shock for a few hundreds of milliseconds after the 
explosion has set in. Moreover, since this energy injection is 
provided by the Vx fluxes, it changes significantly between dif¬ 
ferent progenitors and correlates with the and accretion 
luminosities (at least, during the accretion phase). 


2.4.2. Implementation 

The additional energy deposition, that represents the main 
feature of PUSH, is achieved by introducing a local heating 
term, 2push(^’ (energy per unit mass and time), given by 

X OO 

qUir,E)dE, (4) 


where 


?push(T-E) = 0-0 


1 


Amt, \ tneC^ j Anr^ 


1 



T{r,E), 


(5) 


with 


CTO 


4G^ (nigC^^ 
n (Jicf 


^ 1.759 X lO^'^'^cm^ 


( 6 ) 


being the typical neutrino cross-section, mt, ^ 1.674 x 10“^"^g 
an average baryon mass, and {dLyJdE)l{Anr^) the spectral 
energy flux for any single Vx neutrino species with energy E. 
Note that all four heavy neutrino flavors are treated identically 
by the ASL scheme, and contribute equally to (see the 
factor 4 appearing in Equation (01). 
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The term !F(r, E) in Equation (|5]l defines the spatial location 
where r) is active: 


r(.r,E) 


0 if dsjdr >0 or < 0 

exp(-Ty^(r, £)) otherwise 


(7) 

where Ty^ denotes the (radial) optical depth of the electron 
neutrinos, s is the matter entropy and gy^ y^ the net specific 
energy rate due to electron neutrinos and anti-neutrinos. The 
two criteria above are a crucial ingredient in our description 
of triggering CCSN explosions: PUSH is only active where 
electron-neutrinos are heating (ey^,y^ > 0) and where neutrino- 
driven convection can occur {dsjdr < 0). 

The term Q{t) in Equation (HI determines the temporal be¬ 
haviour of 2push(^’ expression reads 


0 


Q{tj — /^push ^ * 


Un) /^rise 

1 

(Uff ^rise “0 /^rise 

0 


^ — ^on 

Un ^ ^ ton "t" tj-ise 

ton + trise ^ t < toff 
toff ^ t < toff -t" tyfse 
t ^ toff "t tj-fse 


( 8 ) 

and it is sketched in Eigure[3] Note that throughout the article 
we always measure the time relative to bounce, if not noted 
otherwise. 

The cumulative energy deposited by PUSH, Epush, can be 
calculated from the energy deposition rate dEpush/dt as 


EpushCt) — 




df'. 


(9) 


where Ugain is the volume of the gain region. Both these quan¬ 
tities have to be distinguished from the corresponding energy 
and energy rate obtained by IDSA: 




,pdy 


df. (10) 



Figure 3. The function Q{t) determines the temporal behavior of the heating 
due to PUSH. The quantity fon is robustly set by multi-dimensional models. 
We consider a value of 80 ms in our calculations and a value of 120 ms for 
testing, trise ^nd /fpnsh are set by our calibration procedure, spanning a range 
from 50 ms to 250 ms, and from 0 (PUSH off) to ~4, respectively. Since 
we assume that the explosion takes place within the first second after core 
bounce, we use toff = Is. 

The definition of Q{f) introduces a set of (potentially) free 
parameters: 

• ^push is a global multiplication factor that controls di¬ 
rectly the amount of extra heating provided by PUSH. 
The choices of ctq as reference cross-section and of the 


p and T neutrino luminosity as energy reservoir suggest 

^push ~ 1■ 

• ton sets the time at which PUSH starts to act. We re¬ 
late ton to the time when deviations from spherically 
symmetric behavior appear in multi-dimensional mod¬ 
els. Matter convection in the gain region sets in once 
the advection time scale Tadv and the con vective growth 
time scale Tconv satisfy Tadv/Tconv ^ 3 (lEoglizzo et al.l 
120061) . Eor all the models we have explored, this hap¬ 
pens around t = 0.06 - 0.08 s. In the above estimates, 
Tadv = A^shock/A^gain, where Mshock is the accretion rate 
at the shock and Mgain the mass in the gain region, 
and Tconv = fs-v’ where /b-v is the Brunt-Vaisala fre¬ 
quency. Considering that Tconv ~ 4 - 5 ms, we expect 
fon ~ 0.08 - 0.10 s, in agreement with recent multi¬ 
dimensional simulations (see, for exam ple, the ID- 
2D co mparison of the shock position in iBruenn et al.l 

EoBt). 

• frise defines the time scale over which ^(t) in¬ 
creases from zero to kpush- We connect frise with 
the time scale that characterizes the growth of the 
largest multi-dimensional perturbations between the 
shock radius (E.horv) a n d the gain rad i us (Ep ^in) (e.g. 
iJanka & Muelleiiri996h . lEoglizzo et al.l (120061 showed 
that convection in the gain region can be significantly 
stabilized by advection, especially if r^dv/"Tconv only 
marginally exceeds 3, and that the growth rate of the 
fastest growing mode is diminished. Thus, t^se » Tconv 
On the other hand, a lower limit to t^se is represented 
by the overturn time scale, Toveitum. defined as 

^(^ shock " ^gain) /1 1 \ 

Toverturn . v 111 

(V')gain 

where (v)gain is the average fluid velocity inside the gain 
region. In our simulations, we have found Tovertum ~ 
0.05 s around and after ton- In case of a contract¬ 
ing shock, SASIs are al so expected to dev elop around 
0.2 - 0.3 s after bounce (iHanke et al.ll2013h . Hence, we 
assume 0.05 s < frise 5 (0.30 s - fon). 

• foff sets the time after which PUSH starts to be switched 
off. We expect neutrino driven explosions to develop 
for f < 1 s due to the fast decrease of the luminosities 
during the first seconds after core bounce. Hence, we 
fix foff = 1 s. PUSH is not switched off suddenly at 
the onset of the explosion, but rather starts decreasing 
naturally even before 1 s after core bounce due to the 
decreasing neutrino luminosities and due to the rarefac¬ 
tion of the gain region above the PNS. The subsequent 
injection of energy by neutrinos in the accelerating 
shock is qualitatively consistent with multi-dimensional 
simulations, where accretion and explosion can coexist 
during the early stages of the shock expansion. The de¬ 
crease of dEpush/df on a time scale of a few hundreds 
of milliseconds after the explosion has been launched 
makes our results largely independent of the choice of 
foff for explosions happening not too close to it. 

While fon is relatively well constrained and foff is robustly set, 
frise and especially kpush are still undefined. We will discuss 
their impact on th e model and on the explosion properties in 
detail in Sections I3.1.1IITT3] Ultimately, w e wi ll fix them 
using a calibration procure detailed in Section l3/4l 






































Pushing ID CCSNe to explosions: model and SN 1987A 


7 


2.5. Post-processing analysis 

For the analysis of our results we determine several key 
quantities for each simulation. These quantities are obtained 
from a post-processing approach. We distinguish between the 
explosion properties, such as the explosion time, the mass 
cut,or the explosion energy, and the nucleosynthesis yields. 
The former are calculated from the hydrodynamics profiles. 
The latter are obtained from detailed nuclear network calcu¬ 
lations for extrapolated trajectories. 

2.5.1. Accretion rates and explosion properties 

For the accretion process, we distinguish between the ac¬ 
cretion rate at the shock front, Mshock = dM(7?shock)/df, and 
the accretion rate on the PNS, MpNS = dM(/?pNs)/df. In these 
expressions, M{R) is the baryonic mass enclosed in a radius 
R, T^shock is the shock radius, and 7 ?pns is the PNS radius that 
satisfies the conditionp(/?pNs) = 10"gcm“^. 

We consider the explosion time fexpi as the time when the 
shoc k reaches 500 km, me asured with respect to core bounce 
(cf. lUgliano et akl (l2012h '). In all our models, the velocity 
of matter at the shock front has turned positive at that radius 
and the explosion has been irreversibly launched. There is 
no unique definition of 4xni i n the l iterature and some other 
studies (cf. iJanka & Muelie3(ll996l) : lHandv et al.l (120141) 1 use 
different definitions, e.g., the time when the explosion energy 
increases above 10"^* erg. However, we do not expect that 
the different definitions give qualitatively different explosion 
times. 

For the following discussion, we will use the total energy 
of the matter between a given mass shell mo up to the stellar 
surface: 

rma 

FtotaiC^Mo, f) = - I etotai(»J, 0 dm . (12) 

Jm 

M is the enclosed baryonic mass at the surface of the star and 
niQ is a baryonic mass coordinate (0 < mo < M). Ctotai is the 
specific total energy, given by 

^total “ ^int "t" ^kin "t" ^grav 5 ( 13 ) 

i.e., the sum of the (relativistic) internal, kinetic, and gravita¬ 
tional specific energies. For all these quantities we make use 
of the general-relati vistic expressions in the laboratory frame 
(iFischer et al.ll2010h . The integral in Equation (fT^ includes 
both the portion of the star evolved in the hydrodynamical 
simulation and the outer layers, which are considered as sta¬ 
tionary profiles from the progenitor structure. 

The explosion energy emerges from different physical con¬ 
tributi ons (see, for example, the appendix of IScheck et al.l 
(^2006^ and the discussion in lUgliano et al.l (120121) 1. In our 
model, we are taking into account: (i) the total energy of the 
neutrino-heated matter that causes the shock revival; (ii) the 
nuclear energy released by the recombination of nucleons and 
alpha particles into heavy nuclei at the transition to non-NSE 
conditions; (iii) the total energy associated with the neutrino- 
driven wind developing after the explosion up to the end of the 
simulation; (iv) the energy released by the explosive nuclear 
burning in the shock-heated ejecta; and (v) the total (negative) 
energy of the outer stellar layers (also called the “overbur¬ 
den”). We are presently not taking into accout the variation of 
the ejecta energy due to the appearance of late-time fallback. 
This is justified as long as the fallback represents only a small 
fraction of the total ejected mass. 


To compute the explosion energy, we assume that the to¬ 
tal energy of the ejecta with rest-masses subtracted eventually 
converts into kinetic energy of the expanding supernova rem¬ 
nant at f » fexpi- The quantity Ctotai includes the rest mass 
contribution via Cint, see Equation (|2]). Instead, if we want to 
calculate the explosion energy, we have to consider the ther¬ 
mal energy Cth- Therefore, we define the specific explosion 
energy as 

^expl “ ^th + ^kin + ^grav ; (14) 

and the time- and mass-dependent explosion energy for the 
fixed mass domain between mo and M as 

rma 

7/expi(»Jo,0 = - I eexpi(ffi,0dm . (15) 

Jm 

This can be interpreted as the total energy of this region in 
a non-relativistic EOS approach, where rest masses are not 
included. 

The actual explosion energy (still time-dependent) is given 
by 

^expl(0 “ ^expl(^cut(05 0 ■ (^6) 

i.e., for the matter above the mass cut mcut(0- 

To identify the m ass cut, we consider the expression sug¬ 
gested by Bruenn in lEischer et akl (120101) : 

'«cut(0 = m (max(//expi(m, t))j , (17) 

where the maximum is evaluated outside the homologous core 
(m > 0.6 Mq), which has large positive values of the specific 
explosion energy Cexpi once the PNS has formed due to the 
high compression. In the outer stellar envelope, before the 
passage of the shock wave, Cexpi is dominated by the nega¬ 
tive gravitational contribution. However, it is positive in the 
neutrino-heated region and in the shocked region above it. 
Hence, the above definition of mcut locates essentially the tran¬ 
sition from gravitationally unbound to bound layers. The final 
mass cut is obtained for t - funai. 

Our final simulation time ffinai ^ 4.6 s is always much larger 
than the explosion time and, as we will show later, it allows 
Eexpi(0 to saturate. Thus, we consider Eexpi(t - tfinai) as the 
ultimate explosion energy of our models. In the following, if 
we use Eexpi without the time as argument, we mean this final 
explosion energy. 

2.5.2. Nucleosynthesis yields 

To predict the composition of the ejecta, we perform nucle- 
osyn thesis calculations u sing the full nuclear network Win- 
net (IWinteler et akllMT^ . We include isotopes up to ^"Eu 
covering the neutron-deficient as well as the neutron-rich 
side of the valley of ;6-sta b ility. The reaction rates are the 
same as in IWinteler et al.1 (l2012b . They are based on ex¬ 
perimentally known rates where ayailable and predictions 
otherwise. The n-, p-, and al pha-captures are taken from 
iRauscher & ThielemannI (l200(lh . who used known nuclear 
masses where ayailab le and the Finite Range Droplet Model 
(iMoller et al.l Il995bl) for unstable nuclei far from stability. 
The yS-decay rates are from the nuclear database NuDati^ 

We diyide the ejecta into different mass elements of 
10 “^ Mq each and follow the trajectory of each indiyidual 
mass element. As we are mainly interested in the amounts 
of ^®Ni, ^^Ni, ^*Ni, and "^"^Ti, we only consider the 340 in¬ 
nermost mass elements aboye the mass cut, corresponding to 

^ Iittp://www.nndc.bnl.gov/nudat2/ 
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a total mass of 0.34 M©. The contribution of the outer mass 
elements to the production of those nuclei is negligible. 

For t < ffinai, we use the temperature and density evolu¬ 
tion from the hydrodynamical simulations as inputs for our 
network. For each mass element we start the nucleosynthesis 
post-processing when the temperature drops below 10 GK, us¬ 
ing the NSE abundances (determined by the current electron 
fraction Yg) as the initial composition. For mass elements that 
never reach 10 GK we start at the moment of bounce and use 
the abundances from the approximate cr-network at this point 
as the initial composition. Note that for all tracers the fur¬ 
ther evolution of Yg in the nucleosynthesis post-processing is 
determined inside the WtNNET network. 

At the end of the simulations, i.e. t - ffinai, the temperature 
and density of the inner zones are still sufficiently high for nu¬ 
clear reactions to occur (T ^ 1 GK and p ^ 2.5 x lO^g cm“-’). 
Therefore, we extrapolate the radius, density and temperature 
up to fend = 100 s using: 


r(f) — rflnal ^^final 

(18) 

P(0 “ Pfinal ( “ I 

(19) 

\ pfinal / 


T(t) = T[Sf^n^up(t),Ye(t)], 

(20) 


where r is the radial position, v the radial velocity, p the 
density, T the temperature, s the entropy per baryon, and 
Yg the electron fraction of the mass zone. The tempera- 
tu re is calculated at each tim estep using the equation of state 
of iTimmes & Swestvl (l2000h . The prescription in Equations 
(IT8]) -(|201 i corresponds to a free expansion for the density and 
an adiabatic expansio n for the temperature (see, for example, 
iKorobkin et al.l (I2012^ 1. 

3. EITTING AND RESULTS 

To test the PUSH method, we perform a large number of 
runs where we vary the free parameters and explore their im¬ 
pact on the explosion properties. We also analyze in detail 
the basic features of the simulations and of the explosions in 
connection with the properties of the progenitor star. Einally, 
we fit the free parameters in the PUSH method to reproduce 
observed properties of SN 1987A for a progenitor star in the 
range 18-21 Mq. 

3.1. General effects of free parameter variations 

3.1.1. 

The parameter with the most intuitive and strongest impact 
on the explosion is kpush- Its value directly affects the amount 
of extra heating which is provided by PUSH. As expected, 
larger values of kpush (assuming all other parameters to be 
fixed) result in the explosion being more energetic and oc¬ 
curring earlier. In addition, a faster explosion implies a lower 
remnant mass, as there is less time for the accretion to add 
mass to the forming PNS. 

Beyond these general trends with kpush, the detailed behav¬ 
ior depends also on the compactness of the progenitor. Eor all 
16 progenitor models in the 18-21 Mq ZAMS mass range, we 
have explored several PUSH models, varying kpush between 
0.0 and 4.0 in increments of 0.5 but fixing ton = 80 ms and 
Gse = 150 ms. Eor kpush < U none of the models explode 
and for kpush = 1.5 only the lowest compactness models ex¬ 
plode. Eigure |4] shows the explosion energy, the explosion 
time and the (baryonic) remnant mass as function of the pro¬ 
genitor compactness for kpush = 1.5,2.0,3.0,4.0. A distinct 


behavior between low and high compactness models is seen. 
The EC models (^ 1.75 < 0.4) result in slightly weaker and 
faster explosions, with less variability in the explosion energy 
and in the explosion time for different values of kpush- Even for 
relatively large values of kpush, the explosion energies remain 
below 1 Bethe (1 Bethe, abbreviated as 1 B, is equivalent to 
10^' erg). On the other hand, the HC models (^ 1.75 > 0.45) 
explode stronger and later, with a larger variation in the explo¬ 
sion properties. In this case, for high enough values of kpush 
(> 3.0), explosion energies of > 1 Bethe can be obtained. 
The HC models also lead to a larger variability of the remnant 
masses, even though this effect is less pronounced than for the 
explosion time or energy. Eor the values of kpush used here, we 
obtain (baryonic) remnant masses from approximately 1.4 to 
1.9 Mq. The differences of E C and HC models will be inves¬ 
tigated further in Section lJj] 

There are three models with 0.37 < 75 < 0.50 (cor¬ 

responding to ZAMS masses of 18.0 (HC), 18.2 (EC), and 
19.4 M 0 (HC)) which do not follow the general trend. In par¬ 
ticular, we find the threshold value of kpush for successful ex¬ 
plosions to be higher for these models. A common feature of 
these three models is that they have the lowest Ee-core mass of 
all the models in our sample and the highest central densities 
at the onset of collapse. 

The choice of Gse does not affect the observed trends with 
kpush: similar behaviors are also seen for 50 ms < ffise $ 
250 ms. 

3.1.2. ton 

To test the sensitivity of our method to the parameter ton, we 
compute models with k^nsh = 2.0 and ffise = 0.15 s for a very 
large onset parameter, fon = 120 ms. We compare the corre¬ 
sponding results with the ones obtained for ton - 80 ms. As 
expected, the shock revival happens slightly later (with a tem¬ 
poral shift of ~ 30 ms), the explosion energies are smaller (by 
~ 0.05 B) and the remnant masses are marginally larger (by 
0.08 Mq). However, all the qualitative behaviours described 
above, as well as the distinction between high and low com¬ 
pactness models, do not show any dependence on fon- In the 
following, we will always assume fon = 80 ms. 

3.1.3. kpush ^ frise 

In Sections 13.1.11 and 13.1.21 we have investigated the de¬ 
pendency of the model on the single parameters kpush and fon- 
Now, we explore the role of fuse in combination with kpush- Eor 
this, we approximately fix the explosion energy to the canon¬ 
ical value of ~ 1 B for the high compactness models (cor¬ 
responding, for example, to the previously examined models 
with kpush = 3.0 and ffise = 150 ms), and investigate which 
other combinations of kpush and ffise result in the desired ex¬ 
plosion energy. We restrict our explorations to a sub-set of 
progenitor models (18.0 Mq, 18.6 Mq, 19.2 M©, 19.4 M©, 
19.8 Mq, 20.0 Mq, 20.2 Mq and 20.6 Mq) that spans the 
^ 1 . 75 -range of all 16 progenitors. Eigure |5] summarizes the 
explosion energies, explosion times, and remnant masses for 
various combinations of kpush and ffise for progenitors of dif¬ 
ferent compactness. The required constraint can be obtained 
by several combinations of parameters, which lie on a curve 
in the kpush-Gse plane. As a general result, a longer ffise re¬ 
quires a larger kpush to obtain the same explosion energy. This 
can be understood from the different roles of the two parame¬ 
ters: while kpush sets the maximum efficiency at which PUSH 
deposits energy from the reservoir represented by the lu¬ 
minosity, frise sets the time scale over which the mechanism 
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Figure 4. Explosion energies (top), explosion times (middle), and (baryonic) 
remnant mass (bottom) as function of compactness for Arpush 1-5, 2.0, 3.0, 
and 4.0, and fixed frise = 0-15 s for all progenitor models included in this 
studay (ZAMS mass between 18.0 and 21.0 M©). Non-exploding models are 
indicated with Fexpi = —0.5 B in the top panel and are omitted in the other 
panels. 


reaches this maximum. Together, they control the slope of 
0{t) in the rising phase (see Figure|3ll. A model with a longer 
rise time reaches its maximum efficiency later, at which time 
the luminosities have already decreased and a part of the ab¬ 
sorbed energy has been advected on the PNS or re-emitted 
in the form of neutrinos. To compensate for these effects, 
a larger kpush is required for a longer hise. This is seen in 
Figure |6] where we plot the cumulative neutrino contribu¬ 
tion (iipush + Eidsa) and its time derivative for four runs of 


the 18.0 Mq progenitor model, but with different combina¬ 
tions of frise and A:push- Runs with larger parameter values 
require PUSH to deposit more energy (see (iSpush + ^idsa) at 
t ~ fexpi), and the corresponding deposition rates are shifted 
towards later times. Moreover, for increasing values of tnse, 
the explosion time fexpi becomes larger, but the interval be¬ 
tween (ton + frfse) and fexpi decreases. Despite the significant 
variation of ^pu.sh between different runs, the peak values of 
d(£’pu.sh + £idsa)/df at the onset of the shock revival that pre- 
ceeds the explosion are very similar in all cases. 
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Figure 5. Explosion energies (top), explosion times (middle), and (baryonic) 
remnant mass (bottom) as function of compactness for pairs of /cnush and ?rise> 
and for the progenitor models with ZAMS mass 18.0, 18.6, 19.2, 19.4, 19.8, 
20.0, 20.2, and 20.6 M©. 
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Figure 6. Temporal evolution of the total neutrino energy contribution inside 
the gain region (Fpush + ^idsa) (solid lines) and of its time derivative (dashed 
lines), for four runs with the same ZAMS progenitor mass (18.0 Mq), but 
different combinations of PUSH parameters t^sc ipush- For E^rh run, the 
vertical lines correspond to f = ton + trise (long, dashed) and to f^xpl (short, 
dot-dashed). 

3.1.4. toff 

Even though PUSH is active up to Uff + frise ^ 1 s, its en¬ 
ergy deposition reduces progressively on a timescale of a few 
100 ms after the explosion has set in (see Figure |6]l. This 
shows explicitly that the value of toff does not have important 
consequences in our simulations, at least as long as we have 
typical explosion times well below one second. The observed 
decrease of the PUSH energy deposition rat e afte r the launch 
of the explosion will be explained in Section lT3] 

3.2. Contributions to the explosion energy 

In the following, we discuss the contributions to and the 
sources of the explosion energy, i.e., we investigate how the 
explosion energy is generated. This is done in several steps: 
first, we have a closer look at the neutrino energy deposition. 
Then we show how it relates to the increase of the total energy 
of the ejected layers, and finally how this increase of the total 
energy transforms into the explosion energy. For this analysis, 
we have chosen the 19.2 and 20.0 Mq ZAMS mass progeni¬ 
tor models as representatives of the HC and FC samples, re¬ 
spectively. We consider their exploding models obtained with 
ton = 80 ms, frise = 150 ms, and kpush = 3.0. A summary of 
the explosion properties can be found in Table|3] 

The table shows that for both models neutrinos are required 
to deposit a net cumulative energy (Fpush + Fidsa) much larger 
than Fexpi to revive the shock and to lead to an explosion that 
matches the expected energetics. For the two reference runs, 
when the PUSH contribution is switched off (f = foff -i- time), 
the cumulative deposited energy is ~ 4 times larger than Fexpi. 
This can also be inferred from Figure|6]for other runs. That ra¬ 
tio increases further up to ~ 5.5 at f = ffinai, due to the neutrino 
energy deposition happening at the surface of the PNS which 
generates the v-driven wind. According to Equations (|9]l and 
(ITOl i. Fpush and Fidsa are the total energies which are deposited 
in the (time-dependent) gain region. This neutrino energy de¬ 
position increases the internal energy of the matter flowing in 
that region. However, since the advection timescale is much 
shorter than the explosion timescale, a large fraction of this 
energy is advected onto the PNS surface by the accreting mass 
before the explosion sets in, and hence does not contribute to 


Table 3 

Explosion properties for two reference runs 


Quantity 


HC 

LC 

ZAMS 

(Mo) 

19.2 

20.0 

fl.75 

(-) 

0.637 

0.283 

^on 

(ms) 

80 

^rise 

(ms) 

150 

^push 

(-) 

3.0 

^expl 

(ms) 

307 

206 

^remn 

(Mo) 

1.713 

1.469 

^expl (^final) 

(B) 

1.36 

0.57 

^push (^ofF ^rise) 

(B) 

3.51 

1.08 

^idsa (^ofF Mse) 

(B) 

2.76 

1.01 

^idsa (^final) 

(B) 

4.10 

2.11 


Note. — These two runs are used to compare the HC and LC samples. 


the explosion energy. Only the energy deposited by neutrinos 
in the region above the final mass cut will eventually con¬ 
tribute to the explosion energy. 

To identify this relevant neutrino contribution, in Figure [T] 
we show the time evolution of the integrated net neutrino 
energy deposition Ey{ni^ov 0 within the domain above the 
fixed mass = '«cut(ffinai). We choose ni^^i to include all 
the relevant energy contributions to the explosion energy, up 
to the end of our simulations. Despite the significant dif¬ 
ferences in magnitudes, the two models show overall simi¬ 
lar evolutions. If we compare Ey{ni^ov 0 at late times with 
(£'push(foft + frise) + £'idsa(ffinai)) from Tabled we see that it is 
significantly smaller. About two thirds of the energy origi¬ 
nally deposited in the gain region are advected onto the PNS 
and hence do not contribute to the explosion energy. 

In addition to the neutrino energy deposition, in Fig¬ 
ure |7] we also show the variation of the total energy for 
the domain above i.e., AFtotaK^cut’ 0 = ^totaK'Wcut’ 0 - 
Ftotai('«cuo initial), where finitiai is the time when we start our 
simulation from the stage of the progenitor star. The varia¬ 
tion of the total energy can be separated into the net neutrino 
contribution and the mechanical work at the inner boundary, 
AFtotai — Ey + Fmech. We note that in our general relativis¬ 
tic approach the variation of the gravitational mass due to the 
intense neutrino emission from the PNS is consistently taken 
into accout. It is visible in Figure]?] that the net deposition by 
neutrinos makes up the largest part of the change of the total 
energy. The transfer of mechanical energy Fmech is negative 
because of the expansion work performed by the inner bound¬ 
ary during the collapse and the PNS shrinking. However it is 
significantly smaller in magnitude than Ey. 

Next, we investigate the connection between the varia¬ 
tion of the total energy and the explosion energy. In Fig¬ 
ure I?) we show the variation of the explosion energy above 
the fixed mass mlh, i.e. A//expi(»Jcut’ 0 = ??expi(»Jcut> 0 - 
??expi('«cut> initial), together with the relative variation of the 
time-dependent explosion energy, AFgxpKO = £^expi(0 “ 
??expi('«cut> ^initial). It IS obvious from Equations (|2]i, (fT^ . 
and ( fTSi ) that the difference between AHoxpi(mli"i, t) and 
A£totai('w®"(, t) is given by the variation of the integrated rest 
mass energy, AHexpi(«]?;;*(, f) = A£totai(»Jcut> f)-A£mass('«cut> 0- 
In Figure |7] -A£njass(ffJcut’ 0 thus be identified as the 
difference between the long-thin and the short-thick dashed 
lines. We find that the overall rest mass contribution to the 
final explosion energy is positive, but much smaller than the 
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Time post bounce [s] Time post bounce [s] 

Figure 7. Time evolution of the time- and mass-integrated variation of the total energy AEtotai (thin dashed line), of the neutrino net deposition energy Ey (dot- 
dashed line) and of the explosion energy for a fixed domain ATfexpi (thick dashed line), above = mcut(?finai). for the HC (left) and for the LC (right) reference 
nans reported in Table[3] The evolution of the time-dependent explosion energy, AFgxpb is also shown (solid line). Both AT/expi and AFgxpi are computed with 
respect to 7/expi(^^*St’^initial)- The dilference between AFtotaK^^cut’represents the mechanical work, F^ech’ the dilference between AEtotai(^cut’ 
Affexpi(^cut’ 0 represents the released rest-mass energy, — AFmass- 


neutrino contribution. Figure 0 also makes evident the con¬ 
ceptual difference between //expi and Fexpi, and, at the same 
time, shows that //expi(»Jc^t> 0 -> FexpKO for t -> since 
we have chosen = »Jciit(ffinai)- It also reveals that the ex¬ 
plosion energy £expi has practically saturated for f > 1 s, while 
Ey (and, consequently, AFtot and A//expi) increases up to ffinai, 
when is finally ejected. However, this energy provided 
by neutrinos is mostly spent to unbind matter from the PNS 
surface. Thus, the late v-driven wind, which occurs for several 
seconds after 1 s, still increases Fexph but at a relative small, 
decreasing rate. 

To summarize, the variation of the explosion energy above 
can be expressed as 

A//expl('«cut> 0 = AFtotalC'Wcut- 0 “ AFmassC'Wcut- 0 

= t) + £'niech('«cut- 0 - A£'niass('«c;t- 0 • 

( 21 ) 

The quantity -AFniass is positive, but significantly smaller 
than £’v(m|?"j, t). Fmech is negative and also smaller than 
^v('«cut’ 0- Therefore, we conclude that in our models the ex¬ 
plosion energy is mostly generated by the energy deposition 
of neutrinos in the eventually ejected layers, especially within 
the first second after bounce. 

To give further insight, in Figure [8] we show the time evo¬ 
lution of all energies which contribute to the explosion en¬ 
ergy together with the explosion energy itself, for both the 
HC (left panel) and the LC model (right panel). We present 

£'int(»Jcut> 0, 0, £'grav(»Jcut> 0 t), 

which together give a complete decomposition of the explo¬ 
sion energy, i.e., 

f7expl('ncut> 0 = 0 + £'giav(mc;t- 0 

+£.„t(mfi;., t) - £:nass(mcut. 0 ■ (22) 

Compared to Figure |7] we are now not dealing with variations 
any more but with absolute values. Gravitational energy ini¬ 
tially dominates (//expi(fficut’ 0 < 0), meaning that the portion 
of the star above is still gravitationally bound. The HC 


model is initially more bound than the LC model (for exam¬ 
ple, //expl('«cut> f = 0-1 s) ~ -0.54 B, versus //expl('«cm> ^ = 
0.1 s) » -0.40 B, respectively). Before providing positive ex¬ 
plosion energy, neutrinos have to compensate for this initial 
negative binding energy as well as for the negative Lmech- This 
can be seen explicitly by expressing Equation (ISTT i as: 

7fexpl(^cut’ ffinai) ^ hnitial) Ey nal) 

■t"E'inech(ffinai) 5 (23) 

where we have neglected ALmass- 
In the following, we discuss the evolution of the relevant 
ener gies and, in particular, of the rest mass energy (see Sec¬ 
tion 12.21 for the description of the (non-)NSE EOS and of 
the related definitions of the internal, thermal and rest mass 
energies). The innermost part of the ejecta (corresponding 
to ~ 0.15 Mq and ~ 0.07 Mq above for the 19.2 Mq 
and 20.0 Mq model, respectively) is initially composed of in¬ 
termediate mass nuclei (mainly silicon and magnesium). In 
the first part of the evolution, during the gravitational col¬ 
lapse, no significant changes of Eint and Emass are observed 
in Eigure [8] However, when this matter enters the shock, it 
is quickly photodissociated into neutrons, protons, and alpha 
particles. This process increases the rest mass energy, as is 
visible in Eigure [8] between roughly 200 and 300 ms for the 
HC model and between 100 and 200 ms for the LC model. 
At the same time, the release of gravitational energy of the 
still infalling matter and the dissipation of kinetic energy hap¬ 
pening at the shock, together with the large and intense neu¬ 
trino absorption on free nucleons, increase Eint. Later, once 
neutrino heating has halted the collapse and started the ex¬ 
plosion, the expanding shock decreases its temperature and 
free neutrons and protons inside it recombine first into alpha 
particles and then into iron group nuclei. At the same time, 
fresh infalling layers are heated by the shock to temperatures 
above 0.44 MeV, and silicon and magnesium are converted 
into heavier nuclei and alpha particles under NSE conditions, 
leading to an alpha-rich freeze-out from NSE. The production 
of alpha particles, which are less bound than the heavy nu¬ 
clei initially present in the same layers, limits the amount of 
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Time post bounce [s] Time post bounce [s] 

Figure 8. Temporal evolution of the gravitational (thin solid), kinetic (long dashed), negative rest mass (short dashed), internal (dot-dashed), and explosion (solid 
thick) energies above m®"| = mcut(tfinal), for the HC (left) and for the LC (right) reference runs, see Table[3] The internal and the rest mass energy are given with 
respect to the initial rest mass, Fmass.o = Emassif^cLt’ hnitial)- The difference between the internal energy and the rest mass energy represents the thermal energy. 


rest mass energy finally released. Thus, these recombination 
and burning processes liberate in a few hundred milliseconds 
after fexpi an amount of rest mass energy larger but compa¬ 
rable to the energy spent by the shock to photodissociate the 
infalling matter during shock revival and early expansion. We 
have checked in post-processing that the full nucleosynthesis 
network Winnet confirms these results. 

3.3. Explosion dynamics and the role of compactness 

The distributions of the explosion energy and explosion 
time obtained with PUSH, as well as their variations in re¬ 
sponse to changes of the model parameters, suggest a possi¬ 
ble distinction between high and low compactness progeni¬ 
tors. In the following, we investigate how basic properties of 
the models (e.g., the accretion history or the neutrino lumi¬ 
nosities), ultimately connected with the compactness, relate 
to differences in the explosion process and properties. For 
a similar discuss ion in self-consiste nt ID and 2D supernova 
simulations, see iSuwa et ^ (120141) . Again, we choose the 
19.2 and 20.0 ZAMS mass progenitor runs with fuse = 150 ms 
and kpush = 3.0, as representatives of the HC and LC samples, 
respectively. 

In Figure |9] we show the temporal evolution of several 
quantities of interest for both the 19.2 Mq and 20.0 Mq mod¬ 
els, with and without PUSH. The evolution before Ln follows 
the we ll known early shock dyna mics in CCSNe (see, for ex¬ 
ample, iBunow^&^Gosh^ ([19931)). In both models, a few tens 
of milliseconds after core bounce, the expanding shock turns 
into an accretion front, and the mantle between the PNS sur¬ 
face and the shock reaches a quasi-stationary state. In this ac¬ 
cretion phase, Mshock and Mpns are firmly related. However, 
the two different density profiles already affect the evolution 
of the shock. Since P 19 . 2 /P 20.0 ^1-2 outside the shock and 
up to a radius of 2 x 10^ cm (while the infalling velocities 
of the unshocked matter are initially almost identical), Mshock 
(and in turn also MpNs) starts to differ between the two models 
around fpb ~ 30 ms. 

The difference in the accretion rates has a series of imme¬ 
diate consequences. For the HC case, (i) neutrino luminosi¬ 
ties are larger (Figure|9};); (ii) the shock is subject to a larger 
ram pressure (i.e., a larger momentum transfer provided by 
the collectively infalling mass fiowing through the shock). 


and, as visible in the case without PUSH, shock stalling hap¬ 
pens earlier and at a smaller radius (Figure|9l3); (iii) the PNS 
mass grows faster. Since the mass of the PNS at bounce is 
almost identical for the two models (Mpns ~ 0.63 M 0 ), the 
stronger gravitational potential implied by (iii) increases the 
differences in the accretion rates even further by augmenting 
the ratio of the radial velocities inside the gain region (larger 
by 12-15% at f ss Un for the 19.2 Mq case). 

For t > ton, the differences between the two runs amplify as 
a result of the PUSH action. In the LC case, due to the lower 
accretion rate, a relatively small energy deposition by PUSH 
in the gain region (smaller than or comparable to the energy 
deposition by Ve and Vg from IDS A, as visible in Figure flOl i 
is able to revive the shock expansion a few milliseconds after 
ton- Later, the increasing dLpush/dt triggers an explosion in a 
few tens of milliseconds, even before Q{t) reaches its maxi¬ 
mum (Figure |9]2). In the HC case, the energy deposition by 
neutrinos is more intense from the beginning due to the larger 
neutrino luminosities and harder neutrino spectra (Figures [^ 
and|9}l) and due to the higher density inside the gain region. 
However, because of the larger accretion rate, the extra contri¬ 
bution provided by PUSH is initially only able to prevent the 
fast shock contraction observed in the model without PUSH. 
During this shock stalling phase, the accretion rate and the lu¬ 
minosity decrease, but only marginally and very similarly to 
the non-exploding case. When PUSH reaches its maximum 
energy deposition rate (f Un + Lise), the shock revives and 
the explosion sets in (Figure [ 9 ) 2 ). 

In Figure[TT] we plot the ratio of the ram pressure just above 
the shock front (Pram(f^^hock) = calculated at = 

f^shock + 1 km) to the thermal pressure just inside it (Pth(^sbock) 
where R^^ock ~ ^shock - 1 km). In the non-exploding runs (i.e., 
without PUSH), both these pressures decrease with time, but 
their ratio stays always well above unity. On the other hand, 
in runs with PUSH, the more efficient energy deposition by 
neutrinos reduces the decrease of the thermal pressure inside 
the shock. The corresponding drop in the pressure ratio below 
unity determines the onset of the explosion. 

In both runs, once the explosion has been launched, the den¬ 
sity in the gain region decreases and the PUSH energy deposi¬ 
tion rate reduces accordingly. The conversion from an accret- 




























Pushing ID CCSNe to explosions: model and SN 1987A 


13 






Time post bounce [s] Time post bounce [s] 

Figure 9. Temporal evolution of (a) the accretion rate at the PNS and at the shock, (b) the shock, the gain, and the PNS radii, (c) the neutrino luminosities, and 
(d) the neutrino mean energies, for all modeled neutrino flavors. In all panels, we present exploding runs for the 19.2 Mq (red lines) and then 20.0 Mq (blue lines) 
ZAMS mass models obtained with the PUSH parameters reported in Table[3 We also plot the corresponding non-exploding runs obtained by setting kp^sh = 0 
for the 19.2 Mq (light red) and 20.0 Mq (light blue) ZAMS mass progenitor models. 


ing to an expanding shock front decouples Mshock from Mpns ■ 
The latter drops steeply, together with the accretion neutrino 
luminosities (Figures|9h and|9};), while Mshock decreases first 
but then stabilizes around an almost constant (slightly de¬ 
creasing) value. In the case where the shock expansion veloc¬ 
ity is much larger than the infalling matter velocity at /?shock, 
Mshock can be re-expressed as 

Mshock ~ 47r7?^j^Q^j^p(7?shock) Pshock» (24) 

where Vshock = d/?shock/df oc For R > T^shock we have in 

good approximation p{R) oc R^^, and thus 

Mshock ^shock- 

The stationary value of Mshock implies that (5 » 0. Thus, after 
an initial exponential expansion, the shock velocity is almost 
constant during the first second after the explosion. 

Despite the larger difficulties to trigger an explosion, the 
HC model explodes more energetically than the LC model. 
According to the analysis performed in Section 13.21 the dif¬ 
ference in the explosion energy between the HC and the LC 
model depends ultimately on the different amount of energy 
deposited by neutrinos. Since the high compactness model 


requires a larger energy deposition to overcome the ram pres¬ 
sure and the gravitational potential, the total energy of the 
corresponding ejecta (and in turn the explosion energy) will 
be more substantially increased. In addition, after the explo¬ 
sion has been triggered, the larger neutrino luminosities and 
densities that characterize the HC model inject more energy 
in the expanding shock compared with the LC model. 

3.4. Fitting of SN 1987A 

The ultimate goal of core-collapse supernova simulations is 
to reproduce the properties observed in real supernovae. So 
far we have only focused on the dependence of dynamical 
features of the explosion (e.g., the explosion energy) on the 
parameter choices in the PUSH method. However, the ejected 
mass of radioactive nuclides (such as ^®Ni) is an equally im¬ 
portant property of the supernova explosion. Here, we de¬ 
scribe how we calibrate the PUSH method by reproducing the 
explosion energy and mass of Ni ejecta of SN 1987A for a 
progenitor within the expected mass range for this supernova. 

3.4.1. Observational constraints from SN 1987A 

The analysis and the modeling of the observational 
properties of SN 1987A just after the luminosity peak 
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Figure 10. Temporal evolution of the neutrino energy deposition inside 
the gain region from Vg and Vg (dFjdsa/d?, long-thin dashed lines), from 
PUSH (dEpush/d?, short-thick dashed lines), and their sum (solid lines). The 
19.2 M© (H(2) ZAMS mass model is represented in red, while the 20.0 M© 
(LC) ZAMS mass is in blue, with PUSH parameters reported in Table[3] The 
short colored vertical lines show the time of explosion. 



Figure 11. Temporal evolution of the ratio between the ram pressure above 
the shock and the thermal pressure below the shock. The 19.2 M© (HC) 
ZAMS mass model is represented in red, while the 20.0 M© (LC) ZAMS 
mass model is in blue. The PUSH parameters are reported in Table [3] 
Light red and light blue lines represent the corresponding runs without PUSH 
(^push ~ 0)- 


have been the topics of a long s eries of works (e.g. 


WooslevI Il988t lArnett et al.l 11989 

; iShisevama & Nomoto 

19901 Kozma & FranssonI 1998allb 

IBlinnikov et al.l 

200C; 

Fransson & Kozmal 120021 lUtrobin & Chueail 

2005; 


Seitenzahl et alJ 12014 and references therein). They 


provide observational estimates for the explosion energy, 
the progenitor mass, and the ejected masses of ^®Ni, ^^Ni, 
^^Ni, and ^"^Ti, all of which carry rather large uncertainties. 
In Table m the values used for the calibration of the PUSH 
method are summarized. 

The ZAMS progenitor mass is assumed to be between 
18 Mq and 21 M©, corresponding to typical values r eported in 
the literature for the SN 1987A progenitor, see, e.g.. lWooslevl 


(119881) : IShigevama & Nomotol (119901) . For t he explosion en- 
ergy we consider the estimate reported by iBlinnikov et alJ 
(l20()0h . Fexpi = (1-1 ±0.3) X 10^' erg (for a detailed list of 
explosion energy estimates fo r SN 1987A, see for example 
Table 1 in iHandv et al.l (120141) 1. This value was obtained as¬ 
suming ~14.7 Mq of ejecta and an hydrogen-rich envelope of 
~10.3 Mq. The uncertainties in the progenitor properties and 
in the SN distance were taken into account in the error bar. 
The employed values of the total ejecta and of the hydrogen- 
rich envelope are compatible (within a 15% tolerance) with 
a significant fraction of our progenitor candidates, especially 
for Mzams < 19.6 Mq (see Table where the total ejecta 
can be estimated subtracting 1.6 Mq from the mass of the star 
at the onset of the collapse). Explosion models with larger 
ejected mass (i.e., less compatible with our candidate sam- 
ple) tend to have larger e xplosion energies (see, for example, 
lUtrobin & Chugail (120051) 1. Finally, we consider the element 
abund ances for ^®’^^Ni and "^^Ti provided by ISeitenzahl et alJ 
(1201411 . which were obtained from a least squares fit of the de¬ 
cay chains to the bolometric lightcurve. For ^^Ni we use the 
value provided bv iFransson & Kozmal (120021) . 

Table 4 

Observational properties of SN 1987A. 

£expl (l.l±0.3)xl05‘erg 

aiprog 18-21 Mq 

m(5«Ni) (0.071 ± 0.003) Mq 

m(5''Ni) (0.0041 ± 0.0018) Mq 

m(5**Ni) 0.006 Mq 

m(^^Ti) (0.55 ± 0.17) X 10~^ Mq 

Note. — The nucleosynthesis yields are talcen from ISeitenzahl et al.l 
120141) except for ^*Ni which is taken from IFransson & Kozmal 120021) . No 
error estimates were g iven for ^*Ni. The explosion energy is adapted from 
IBlinnikov et al.l 120001) . For the progenitor range we chose typical val- 
ues found in the literature, see e.g. IShigevama & Nomoto l 11990 1): IWooslevI 

CMa. 


3.4.2. Fitting procedure 

We calibrate the PUSH method by finding a combination 
of progenitor mass, kpush, and fuse which provides the best 
fit to the all observational quantities of SN 1987A mentioned 
above. The weight given to each quantity is related to the un¬ 
certainty. For example, due to the large uncertainty in the '^"^Ti 
mass, this does not provide a strong constraint on selecting 
the best fit. 

Figure [12] shows the explosion energy and ejected mass of 
^®Ni, ^^Ni, ^*Ni, and "^"^Ti for different cases of kpush and frise 
and for four select HC progenitors used to calibrate the PUSH 
method. We do not consider the FC progenitors, because of 
their generally lower explosion energies, see Figure ID The 
different cases of kpush and frise span a wide range of explo¬ 
sion energies around 1 Bethe. For all parameter combinations 
shown, at least one progenitor in the 18-21 Mq range fulfills 
the requirement of an explosion energy between 0.8 Bethe and 
1.4 Bethe. There is a roughly linear correlation between the 
explosion energy and the synthesized ^®Ni-mass. However, 
this correlation is not directly compatible with the observa¬ 
tions, as the ejected ^®Ni is systematically larger than expected 
(up to a factor of ~ 2 for models with an explosion energy 
around 1 Bethe). There is a weak trend that models with 
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Eexpi[10®^ erg] Egxpi [10®^ erg] 



Figure 12. Ejected mass of (top left), (top right), (bottom left), and ^Ti (bottom right) and explosion energy for four representati ve HC progenitor 
models. Five combinations of ipush and fuse are shown, each with a different symbol. The error bar box represents the observational values from ISeitenzahl et all 
<201411 (for 56.57jyfj gpij 44^1^ 2 nd from IFransson & Kozmal <20021) (for ^*Ni). No error bars are reported for ^*Ni. 


higher f^se tend to give lower nickel masses for given explo¬ 
sion energy. Among the parameter combinations that produce 
robustly high explosion energies (i.e., kpush > 3), kpush = 3.5 
with the high value of tnse of 200 ms gives the lowest ^®Ni 
mass for similar explosion energies, but still much too high. 

Our simulations can be reconciled with the observations by 
taking into account fallback from the initially unbound mat¬ 
ter. Since we do not model the explosion long enough to see 
the development of the reverse shock and the appearance of 
the related fallback when the shock reaches the hydrogen-rich 
envelope, we have to impose it, removing some matter from 
the innermost ejecttQ. With a value of ~ 0.1 Mq we can match 
both the expected explosion energy and ^®Ni ejecta mass, see 
Figure [13] In this way we have fixed the hnal mass cut by 
observations. However, we point out that we are able to iden¬ 
tify the amount of late-time fallback only because we also 
have the dynamical mass cut from our hydrodynamical sim¬ 
ulations. This is not possible in other methods such as pis¬ 
tons or thermal bombs. Our value of ~ 0.1 Mq of fallback in 

'll Note that we did not modify the explosion energy due to the fallback. 
This is based on the expectation that at the late time when fallback forms, the 
explosion energy is approximately equally distributed among the total ejected 
mass, which is about two orders of magnitude higher than our fallback mass. 


SN 1987A will b e fur ther discussed and compared with other 
works in Section l43] 

The observed yield of ^®Ni provides a strong constraint on 
which parameter combination would fit the data. From the ob¬ 
served yields of ^^Ni and ^^Ni, only the 18.0 and 19.4 progen¬ 
itors remain viable candidates. Without fallback our predicted 
'^‘*Ti yields are compatible with the observed yields (see Fig¬ 
ure Ell. However, if we include fallback (which is needed to 
explain the observed Ni yields), ^'^Ti becomes underproduced 
compared to the oberved value. Since this behavior is true 
for all out models, we exclude the constraint given by '^Ti 
from our calibration procedure. From the considered param¬ 
eter combinations, we obtained the best ht to SN 1987A for 
the 18.0 Mq progenitor model with kpush = 3.5, fuse = 200 ms, 
and a fallback of 0.1 Mq. These parameters are summarized 
in Table |5] In Figure [T4| we show the temporal evolution of 
the accretion rates, of the relevant radii, and of the neutrino lu¬ 
minosities and mean energies for our best fit model. For com¬ 
parison purposes, we present also the results obtained for the 
same model without PUSH. Note that in this non-exploding 
case the Ve and Ve luminosities stay almost constant over sev¬ 
eral ~ 100 ms after core bounce, despite the decreasing accre¬ 
tion rate. This is due to the relatively slow variation of Mpns 
( for example, compared with the variation obtained in the 19.2 
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Figure 13. Same as Figure fTTl but assuming 0.1 Mq fallback. Note the different scale for ^®Ni and ^*Ni compared to Figure fTTl 
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• 19.2 
19.4 

• 20.6 
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Mq model. Figure |9ll and due to the simultaneous increase of 
the PNS gravitatio nal potential, propo rtional to Mpns/^pns 
( see, for example. iFischer et al.ll200^ . A summary of the 
most important results of the simulations using this param¬ 
eter set for the different progenitors in the 18-21 Mq window 
is given in Table |6] For the remnant mass and for the ^®Ni 
yields of our best-fit model, we provide both the values ob¬ 
tained with and without assuming a fallback of 0.1 Mq. 

Table 5 

Parameter values for best fit to SN 1987A. 


^push 

^rise 

^on 

^ofF 

(-) 

(ms) 

(ms) 

(s) 

3.5 

200 

80 

1 


Note. — We identified the 18.0 M© model as the progenitor which fits 
best, whereas we had to impose a late-time fallback of 0.1 M©. 


3.5. M and Ti yields, progenitor dependence 

Figures[T2]and[T3]show that the composition of the ejecta is 
highly dependent on the progenitor model, especially for the 
amount of ^^Ni and ^^Ni ejected. From the four HC progeni¬ 
tors shown, two (18.0 Mq and 19.4 Mq) produce a fairly high 


Table 6 

Summary of simulations for /cpush = 3.5 and = 200 ms. 


ZAMS 

^expl 

^expl 

•‘remnant 

remnant 

M(^*Ni) 

(Mo) 

(Bethe) 

(s) 

(Mo) 

(Mo) 

(Mo) 

18.0 

1.092 

0.304 

1.563 

1.416 

0.158 

18.2 

0.808 

0.249 

1.509 

1.371 

0.110 

18.4 

1.358 

0.318 

1.728 

1.549 

0.144 

18.6 

0.702 

0.239 

1.529 

1.388 

0.090 

18.8 

0.721 

0.236 

1.522 

1.382 

0.093 

19.0 

1.366 

0.317 

1.716 

1.54 

0.161 

19.2 

1.356 

0.318 

1.724 

1.546 

0.152 

19.4 

1.15 

0.326 

1.608 

1.452 

0.158 

19.6 

0.371 

0.230 

1.584 

1.433 

0.04 

19.8 

0.661 

0.225 

1.523 

1.383 

0.088 

20.0 

0.613 

0.222 

1.474 

1.342 

0.085 

20.2 

0.379 

0.224 

1.554 

1.408 

0.039 

20.4 

0.743 

0.263 

1.674 

1.506 

0.094 

20.6 

1.005 

0.277 

1.781 

1.592 

0.141 

20.8 

0.959 

0.277 

1.764 

1.578 

0.135 

21.0 

1.457 

0.316 

1.733 

1.554 

0.198 

18.0 (lb) 

1.092 

0.304 

1.663 

1.497 

0.073 

Note. — For 

the model 

18.0 (fb), which is 

our best fit 

to SN 1987A, 


have included 0.1 M© of fallback, determined from obervational constraints. 
See the text for more details. 
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Time post bounce [s] 




Time post bounce [s] 


Figure 14. Same as in Figure|9] but for the SN 1987A best fit model: 18.0 Mq progenitor, with = 80 ms, tnse = 200 ms, and <:push = 3.5. 


amount of those isotopes, while the other two (19.2 Mq and 
20.6 Mq) do not reach the amount observed in SN 1987A. A 
thorough investigation of the composition profile of the ejecta 
reveals that ^^Ni and ^^Ni are mainly produced in the slightly 
neutron-rich layers (Yg < 0.5), where the alpha-rich freeze- 
out leads to nuclei only one or two neutron units away from 
the N - Z line. A comparison of the and composition pro¬ 
files for the 18.0 Mq and the 20.6 Mq progenitors is shown in 
Figure[T5] For the 18.0 Mq model, the cutoff mass is 1.56 Mq 
and a large part of the silicon shell is ejected. In this shell, the 
initial matter composition is slightly neutron-rich (due to a 
small contribution from ^®Fe) with Y^ — 0.498 (dotted line in 
top left graph) and the conditions for the production of ^^Ni 
and ^^Ni are favorable. The increase in Y^ around 1.9 Mq 
marks the transition to the oxygen shell. The same transition 
for the 20.6 Mq model happens around 1.74 Mq, i.e., inside 
the mass cut. Therefore, this mod el ejects less ^^Ni and ^*Ni 
(see also iThielemann et al.lll99()h . In all our models, ^"^Ti is 
produced within the innermost 0.15 Mq of the ejecta (see Fig- 
urefTSll. Since we assume 0.1 Mq fallback onto the PNS, most 
of the synthesized ^"^Ti is not ejected in our simulations. 


4. IMPLICATIONS AND DISCUSSION 
4.1. Sensitivities of nucleosynthesis yields 


While post-processing the ejecta trajectories for nucleosyn¬ 
thesis, Ye is evolved by the nuclear network independently of 
the hydrodynamical evolution. This leads to a discrepancy at 
later times between the electron fraction in the initial trajec¬ 
tory (Tg^‘*™) and in the network (T™‘^). In order to estimate 
the possible error in our nucleosynthesis calculations arising 
from this discrepancy, we have performed reference calcula¬ 
tions using Ye^'^'°(t - ffinai) instead of = 10 GK) as 

a starting value for the network (see Section |2. 5. 21 l. The re¬ 
sults are shown in Figure fTSl for two progenitors: 18.0 Mq 
and 20.6 Mq. The label “standard” refers to the regular case 
which uses Y'f^'^^°{T - 10 GK) as input. The calculation us¬ 
ing Te^*°(f = ffinai) as input is labeled “alternative” and is 
represented by the dashed lines. The point in time at which 
the Ye prohle is shown is indicated by the supplements “in¬ 
put” (before the first timestep) and “hnal” (at t - 100 s). The 
corresponding nuclear compositions of the ejecta, each at the 
hnal calculation time of 100 s, are shown in the bottom pan¬ 
els. For the alternative Ye prohle of the 18.0 Mq progenitor 
(top left) the minimum around 1.59 Mq disappears, leading to 
an increase in ^®Ni in this region at the expense of ^’Ni and 
^^Ni (bottom left). For the 20.6 Mq progenitor the situation 
is similar, with only a very small region just above 1.8 Mq 
showing signihcant differences. In general, we observe that 
the uncertainties in Ye in our calculations are only present up 
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Figure 15. Electron fraction profiles (top) and nuclear compositions at 100 s (bottom) above the mass cut for the 18.0 Mq (left) and the 20.6 Mq (right) 
progenitors with the parameters /ipush = 3.5 and fuse = 200 ms. The electron fraction is plotted for two different times in the network: the input values for the first 

timestep (“inpuf’) and the value after post-processing (“final”). The dashed lines in all panels correspond to the alternative case, where = 4.6 s) is taken 

as the initial electron fraction in the network, whereas the solid lines represent the standard case (using = 10 GK)). 


to 0.05 Mq above the mass cut. The resulting uncertainties in 
the composition of the ejecta are very small or even inexistent 
in the scenarios where we consider fallback. 

The radioactive isotope "^'^Ti can be detected in super¬ 
novae and supernova remnants. Several grou ps have used 
differ e nt techniques to estimate th e ' ^"^Ti yield (IChugai et al.l 
1 9971; iFransson & Kozm^ |2002|; iJerkstrand et al.j 201 U 


Larssraet^ 201_1^ &ebenevet al.ll2012t iGrefenstette et akl 


2014t ISeitenzahl et al.l l2014h . The inferred values span a 
broad range, (0.5 - 4) x 10“^ Mq . Traditional supernova 
nucleosynthesis ca l culatio ns (e.g. iThielemann et al.l Il996t 
IWooslev & Weav^ 119951) typically predict too low '^Ti 
yields. Only ve r y few models predict high '^‘^Ti yields: 
IThielemann et al.l (il990h report '^‘*Ti yields around 10 
and above in t he best fits of t heir a rtificial SN explosions 
to SN 1987A. iRauscher et al.l (I2002^ argue that the yields 
of ^®Ni and ^"^Ti are very sensitive to the “final mass cut” 
(as we have shown, too), which is often determined by 
fallback. Ejecta in a supernova may be subject to convective 
overturn. To account for this, we can assume homogeneous 
mixing in the inner layers up to the outer boundary of the 
silicon sh ell before cutting off the f allback material (see, for 
example, lUmeda & Nomotc^ (120021) and references therein). 
For our best-fit model, the ejected "^^Ti mass increases to 
2.70 X 10“^ Mq, if this prescription is applied. Comparing 


to the previous yield of 1.04 x 10^^ Mq, we observe that 
the effect of homogeneous mixing is considerable, but not 
sufficient to match the observational values. The ejected 
56-58]sjj masses also show a slight increase. However, there 
are also uncertainties in the nuclear physics connected to the 
production and destruction of "'^Ti. The final amount of pro¬ 
duced "^^Ti depends mainly on two reactions: "^°Ca(Q', y)"^"^Ti 
and py'^V. Recent measurements of the '*"^Ti(Q', p)"^^V 

reaction rate within the Gamow window concluded that 
it may be considerably smaller than preyious theoretical 
predictions (iMargerin et al.l [^14ft . In this study, an upper 
limit cross section is reported that is a factor of 2.2 smaller 
than the cross section we haye used in our calculations 
(at a confidence leyel of 68%). Using this smaller cross 
section for the "^"^Ti(a;, p)"^^V reaction, our yield of ejected 
"^"^Ti for our best-fit model (18.0 Mq progenitor, kpush = 3.5, 
(rise = 200 ms) rises to 1.49 x 10“^ Mq with fallback and 
5.65 X 10“^ Mq without fallback. This corresponds to 
a relatiye increase of 43% with fallback and 48% without 
fallback. If we include both the new cross section and 
homogeneous mixing, the amount of ^"^Ti in the ejecta is 
3.99 X 10“^ Mq including fallback. This yalue, howeyer, is 
still below the expected yalue deriyed from obseryations, but 
within the error box. 
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4.2. Wind ejecta 

In the analysis of the nucleosynthesis yields above, we have 
used a mass resolution of 0.001 M© for the tracers. This is 
too coarse to resolve the ejecta of the late neutrino-driven 
wind. Note that in our best-fit approach, where no mixing 
is assumed, none of the neutrino-driven wind is ejected be¬ 
cause it is part of the fallback. Nevertheless, in the follow¬ 
ing we report briefly on the properties of the wind obtained 
by our detailed neutrino-transport scheme. For our best-fit 
model, the 18.0 Mq progenitor, at ffinai we And an electron 
fraction around 0.32, entropies up to 80 kg per baryon, and 
fast expansion velocities (~ 10^ cm/s). Similar conditions are 
also found for the other progenitors. They are not sufficie nt 
for a full r-process (see, for example, iFarouqi et al.l (l2010h l. 
On the other hand, we have found that the entropy is still in¬ 
creasing and the electron fraction still decreasing in the fur¬ 
ther evolution. The high asymmetries are only obtained if we 
include the nucleon mean-fiel d interaction potentials in th e 
neutrino charged-current rates (iMartlnez-Pinedo et al.ll201^ . 
However, they are much higher than found in ot her long-term 
simulations which also include these pote ntials (iRoberts et alJ 
120121 IMartlnez-Pinedo et al.ll201^ l2014h . This could be re¬ 
lated to the missing neutrino-electron scattering in our neu¬ 
trino transport, which is an important source of thermalization 
and down-scattering, especial ly for the high energ y electron 
antineutrinos at late times, see lFischer et al.l (120121) . More de¬ 
tailed comparisons are required to identify the origin of the 
found differences which will be addressed in a future study. 

4.3. Amount of fallback 

To reconcile our models with the nucleosynthesis obsery- 
ables of SN 1987A we need to inyoke 0.1 M© of fallback (see 
Section [3.4.2b . The yariation in the amount of synthesized Ni 
isotopes between runs obtained with different PUSH parame¬ 
ters (FigurefT2li suggests that a smaller t^se (and consequently 
smaller kpush) could also be compatible with SN 1987A ob- 
seryables, if a larger fallback is assumed. On the one hand, 
assuming that t^se ranges between 50 ms and 250 ms, fall¬ 
back for the 18.0 M© model compatible with obseryations 
is between 0.14 M© (for frise = 50 ms) and 0.09 M© (for 
Irise = 250 ms). On the other hand, if the amount of fall¬ 
back has been fixed, the obseryed yields (especially of ^®Ni) 
reduce the uncertainty in frise to < 50 ms. 

Our ch oice of 0.1 M© is co mpatible with the fallback ob¬ 
tained by lUgliano et alJ (120121) in exploding spherically sym¬ 
metric models for p rogenitor s t ars in the same ZAMS mass 
window. Moreoyer, ICheyalieJ (Il989h estimated a total fall¬ 
back around 0.1 M© for SN 1987A, which is supposed to be 
an unusually high yalue compared to “normal” type II super- 
npyae. Recent mu lt i-dimensional num erical simulations by 
iBernal et alJ (120131) : iFraiia et ^ (1201 4ft confirmed this sce¬ 
nario and furthermore showed that such a hypercritical accre¬ 
tion can lead to a submergence of the magnetic held, giying 
a natural explanation why the neutron star (possibly) born in 
SN 1987A has not been found yet. 

4.4. Compact Remnant of SN 1987A 

From the obseryational side, the compact remnant in 
SN 1987A is still o b scure. From the neut rino signal (see , e.g., 

I Arnett et ^ (ll989ft : lKoshibal (1 1992ft and IVissaM (l2()15ft for a 
recent detailed analysis) one can conclude that a PNS star was 
formed and that it lasted at least for about 12 s. The mass cut 
in our calibration run is located at an enclosed baryon mass of 


1.56 M© without fallback. If we include the 0.1 M© of late¬ 
time fallback required to fit the obseryed nickel yields and the 
explosion energy, we haye a final baryonic mass of 1.66 M©. 
For the employed HS(DD2) EOS this corresponds to a grayi- 
tational mass of a cold neutron star of 1.42 M© (without fall¬ 
back) or 1.50 M© (wit h fallback). The CCSN sim ulations with 
artificial explosions of iThielemann et aTI (11990ft . where a Anal 
kinetic energy of 1 Bethe was obtained by hand and where the 
mass-cut was deduced from a ^®Ni yield of (0.07 + 0.01) M©, 
lead to a similar baryonic mass of (1.6 + 0.045)Mo. These au¬ 
thors also wrote that uncertainties in the stellar models could 
increase this yalue to 1.7 M© which would also be fully com¬ 
patible with our result. 

The prediction of the neutron star mass has important con¬ 
seq uences. From the obsery ations of iDemorest et al.l (12010ft 
and lAntoniadis et alJ (12013ft it follows that the maximum 
grayitational mass of neutron stars has to be aboye two so¬ 
lar masses. The maximum mass of the HS(DD2) EOS is 
2.42 M©, corresponding to a baryonic mass of 2.92 M©. If 
the compact remnant in SN 1987A was a black hole, and not 
a neutron star, it means that at least ~ 0.5 M© of additional ac¬ 
creted mass were required, if we just take the two solar mass 
limit. If we use the maximum baryonic mass of HS(DD2) we 
eyen haye to accrete ~1.3 M© of additional material. Obyi- 
ously, if such a huge amount of material would be accreted 
onto the neutron star, our predictions for the explosion energy 
and the nucleosynthesis would not apply any more. 

Neyertheless, we haye the impression that it would be diffi¬ 
cult to fit the SN 1987A obseryables and obtain a black hole as 
the compact remnant at the same time. Eor spherical fallback, 
it is certainly excluded. The only possibility could be a highly 
anisotropic explosion and aspherical accretion, which we can¬ 
not address with our study. To show if such a scenario can be 
realized remains a task for future multi-dimensi onal studies. 
In the 2D simulations of lYamamoto et all (12013ft the remnant 
mass is decreasing with the explosion energy and an explo¬ 
sion energy aboye 1 Bethe would res ult in neutron stars be- 
low ~ 2 M© baryonic mass. Note that iKifonidis et af] (120061) 
already came to the same conclusion that the formation of a 
black hole in SN 1987A “is quite unlikely”, based on 2D sim¬ 
ulations with a 15 M© progenitor. 

Another possibility was proposed by IChan et ^ (12009ft . 
These authors argued that the time delay of ~ 5 s obseryed 
for the neutrino signal by the IMB detector could be related to 
a collapse to a quark star. Due to the proposed faster neutrino 
cooling of quark stars, this would giye a natural explanation 
why it has not been obseryed until today. The end of our simu¬ 
lations is also around 5 s, thus we can make statements about 
the conditions at which the phase transiti on to quark matter 
took place in SN 1987A, if the scenario of IChan et al.l (12009ft 
was true. We haye a central mass density of 4.56 x lO'"^ g/cm^ 
corresponding to = 0.272 fm“^ or = 1.83 a tem¬ 
perature of 23.2 MeV, and an electron fraction of 0.24. Some 
simplified models for quark matter predict that the phase tran¬ 
sition in symmetric matter is shif ted to higher densiti es com¬ 
pared with supernoya conditions (lEischer et alJl201 H) . Under 
that hypothesis, a phase transition around 2 po and 20 MeV 
cannot be excluded. 

A simpler explanation is giyen by the possibility that a pul¬ 
sar in the SN 1987A remnant is simply not (yet) obseryable. 
iQgelman & Alpatl (12004ft : iGrayes et al.l (12005ft showed that 
the non-detection of any compact remnant puts important lim¬ 
its on the magnetic held the NS can haye (either unusually low 
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or very high, in the real m of magnetars). F urthermore, for 
both cases (NS and BH) [Graves et al.l (120051) put severe con¬ 
straints on currently ongoing accret ion scenarios , e.g., s pher- 
ical accretion is almost ruled out. iGraves et al.l (l2005h con¬ 
clude that “it seems unlikely that the remnant of SN 1987A 
currently harbors a pulsar”. Our simulations would be in line 
with the option of a neutron star with a very low magnetic 
field or with a “normal” magnetic field which is still (partly) 
buried in the crust due to the late time fallback, similar to what 
is observed for neutron stars in binary systems. In this respect, 
recent high-resolution radio observations of the remnant indi- 
cate the presence of a comp act source or a pulsar wind nebula 
(IZanardo et al.ll20131120141) . Future observations will be able 
to clarify the nature of this emission. 

4.5. Correlations 

As a byproduct of exploring the 18-21 Mq window and the 
fitting procedure to SN 1987A we have found interesting cor¬ 
relations between different quantities, which we will discuss 
here. In Figure[T 6 ] we plot the explosion energy, the explosion 
time, and the (baryonic) remnant mass as function of the pro¬ 
genitor compactness. The results obtained with the calibrated 
runs indicate a general trend with progenitor compactness for 
Fexpi- The explosion time, fexpU is almost constant within each 
the LC and the HC group, while the difference between the 
two groups is related to the difference be twee n how LC and 
HC models explode (discussed in Section ISJl l. The remnant 
mass increases with compactness, as expected. Nevertheless, 
we notice significant deviations from the described trends; for 
Eexpi and fexpi in the HC sample, for Mrem mainly in the LC 
sample. 

Figure[T7]shows explosion times and explosion energies for 
all the exploding runs in our sample. We can identify a corre¬ 
lation between fexpi and Lexpi for a given progenitor: the larger 
fexpi the lower is Fexpi- This correlation is more pronounced 
for the HC models than for the LC models. It means that the 
explosion in PUSH cannot set in too late, if the observed ex¬ 
plosion energy should be achieved. 


4.6. Heating efficiency and residence time 

In the context of CCSNe, the heating efficiency 77 is often 
defined as the ratio between the volume-integrated, net energy 
deposition inside the gain region and the sum of the Vg and Vg 
luminosities at infinity: 



Fvg "f" Ey^ 


(26) 


see, e.g., iMurphv & Burrow^ (120081):iMarek & Jankal (l2009h: 

iMuller et alJ(l2012ah : ICouch & O’Connoii(l20l4 :IS uwa et al] 

(120131) . ill non-exploding, spherically symmetric simula¬ 
tions, 77 usually rises within a few tens of milliseconds af¬ 
ter core bounce and reaches its maximum around 77 ~ 0.1 at 
f ^ 100 ms, when the shock approaches its maximum radial 
extension. As soon as the shock starts to recede and the vol¬ 
ume of the gain region decreases, 77 diminishes quickly to a 
few percents (see, for example, the long-dashed lines in Fig¬ 
ure [T^. 

In multi-dimensional simulations, where the shock contrac¬ 
tion is delayed or even not happening, energy deposition is 
expected to be slightly more efficient (77 ~ 0.10 - 0.15 at max¬ 
imum) and to decrease more slowly, within a few hundreds of 
milliseconds after bounce or at the onset of an explosion (see. 


for ex ample, |Mur£h^_&^unow^j200S iMuller et alJ 12012at 
ICouchl 120131 ICouch & O’Connoill2014ft . These differences 
arise not only because the gain region does not contract, but 
also because neutrino-driven convection efficiently mixes low 
and high entropy matter between the neutrino cooling and the 
heating regions below the shock front. Furthermore, convec¬ 
tive motion and SASIs are expected to increase significantly 
the residence time of fluid particles inside the gain region dur- 
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Figure 16. Explosion energies (top), explosion times (middle), and (bary¬ 
onic) remnant mass (bottom) as function of compactness for the PUSH pa¬ 
rameters of our best-fit model (/cpush = 3.3 and ?rjse = 0.15 s) for all progen¬ 
itors in the 18-21 M© window. HC models are denoted by a red cross, LC 
models by a blue plus. Our best-fit model for SN 1987A is highlighted by a 
filled triangle. 
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Figure 17. Explosion energy Eexpi versus explosion time texpi for Iho 
progenitors in the 18-21 M© range and for different combinations of /cpush 
and trise, however only the exploding models are included. HC models ai‘e 
indicated by a triangle, LC models by a circle. The best fit model is indicated 
by a cross. The different colors distinguish different progenitors. 


ties are uni vocally defined by the cooling of the PNS and 
by the accretion rate history, we increase the energy gain by 
acting on the neutrino heating efficiency. This effect can be 
made visible by defining a heating efficiency that takes the 
PUSH contribution into account, ? 7 tot: 


^tot ~ n ^push ~ 





(27) 


In Figure [181 we plot rfiot as a function of time for our 
SN 1987A calibration model, with PUSH (kpush = 3.5) and 
without it (kpush = 0). We first notice that the heating effi¬ 
ciency provided by Vg and Vg can differ between exploding 
(short-thick dashed lines) and non-exploding models (long- 
thin dashed lines). In the case of the exploding model, PUSH 
provides an increasing contribution to T/tot. It continues to in¬ 
crease steeply up to f a; ton + frfse. but also later, up to f a; fexpi, 
due to the shock expansion preceding the explosion. Thus, 
the increasing heating efficiency in our spherically symmet¬ 
ric models can be interpreted as an effective way to include 
average residence times longer than the advection timescale. 






Figure 18. Neutrino heating efficiency for the SN 1987A best fit model: 
18.0 M© model with ?on = 80 ms, fj-ise = 200 ms, and Aipush = 3.5. The 
solid lines represent the total efficiency (i.e., due to Vg and Vg absorption and 
due to PUSH), the short-thick dashed lines the efficiency only due to Vg and 
Vg absorption. For comparison, the heating efficiency of the corresponding 
non-exploding model (/cpush = 0 ) is also presented (long-thin dashed lines). 


Figure 19. Average and maximum heating efficiencies, calculated between 
t = ton and t = ^expi for the runs obtained with the fitted parameters, Table[5] 
and plotted as a function of the progenitor compactness ^ 1 . 75 . The black 
crosses and the red triangles refer to the average and the maximum efficiency 
due to Vg and Vg {rj), while the blue stars and the magenta squares to the 
average and the maximum total efficiency (//tot), including also the PUSH 
contribution. 


ing which they are subject to i ntense neutrino he ating (see, 
e.g., iMurphv & Burrov^l2008UHandv et al.ll2014^ . Since the 
increase of the particle internal energy is given by the time 
integral of the energy absorption rate ove r the residence time , 
this translates to a larger energy variation (iHandv et al.l2014l) . 

In spherically symmetric models, the imposed radial 
motion does not allow the increase of the residence time. 
This constraint limits the energy gain of a mass element 
traveling through the gain region. In models exploded 
using the light-bulb approximation, a large enough internal 
energy variation is provided by increasing the neutrino 
luminosity above a critical value, which depends on the 
mass accretion rate and on the di mensionality of the 


mode l (e.g.. iBurrows & GoshvI 19931 lYamasaki & Yamadal 


iwakarni et al.l 120081: Murphy & Burrows ~ |20()8t' 
Iwakami et al.l ^0^ iNor dhaus et al.l 1201 o[~ Hankeetal.l 


2001 


20121 ICouchlf2013l: iDolence et all 120131: IHandv et all 12014 

Suwa et al.ll2014l) .~Since in our model the neutrino luminosi- 


In Figure[T9] we collect the average and the maximum heat¬ 
ing efficiencies, for all the models obtained with the set of 
parameters resulting from the fit procedure (Table [S]). Both 
the average and the maximum values are computed within the 
interval ton < t < fexpi. We plot them as a function of the com¬ 
pactness and we distinguish between tj and T/tot. The maxi¬ 
mum of 7 / is usually realized att ^ ton, while the maximum of 
7/tot is reached around t « fexpi (see also Figure [TSll. Since the 
explosion sets in later for HC models, when fexpi > fon + frfse. 
the PUSH factor 0 has time to rise to kpush for these models. 
This increases not only the maximum but also the average T/^t 
compared with the LC cases. We notice that all four quan¬ 
tities show a correlation with fi. 75 , but much weaker in the 
case of 7 / than in the case of 7/tot. Moreover, in the HC region, 
we recognise deviations from monotonic behaviors which re¬ 
produce the irregularities already observed in the explosion 
properties. 
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4.7. Alternative measures of the explosion energy 

In the following, we discuss alternative measures of the 
explosion energy used in the literature for reasons of com¬ 
parison. We investigate their behaviors at early simulation 
times and their general ra te of convergence. The diagnos¬ 
tic energy E^(t), see e.g. iBruenn et alJ (120131) . is given by 
the integral of the specific explosion energy Cexpi over re¬ 
gions wher e it is positive (again, excluding the PNS core, 
see Section l2.5.1b . The quantity E^{t) is often used in multi¬ 
dimensional simulations as an estimate of the explo s ion en ¬ 
ergy at ear l y sirn u lation times, see e.g . iBuras et al.l (120061): 
ISuwa et al.l(f2010l):lMulle r et al.l(l^l2bl) : ICouch & O’Connoii 
(l20l4 : lTakiwaki et al.l(l20T^ 


The overburden fiovCO. see lBruenn et al.l(l2013h . is given by 
the integral of the specific explosion energy of the still grav¬ 
itationally bound regions between the expanding shock front 
and the surface of the progenitor star. If we define E^^{t) as 
the sum of the overburden and of the diagnostic energy, we 
recover a measure of the explosion energy equivalent to the 
one defined in Equation ( fTSl i: 


£expl(0 = C(0 = E^(t) + £ov(f)- 


(28) 


For long enough simulation times, all matter above the mass- 
cut should get positive specific explosion energies, and thus 
the overburden should approach zero and the diagnostic en¬ 
ergy should become equal to the explosion energy £expi(0- 
Finally, an upper limit for the explosion energy is obtained 
by also takin g into account the “residual recombination en¬ 
ergy” £'rec(0 ( Bruenn et al.ll2013h : 


Eov.r(t) = 'E'ov(0 + £'iec(0 , 


(29) 


where Free)?) is the energy that would be released if all 
neutron-proton pairs and all "^He recombined to ^®Ni in the 
regions of positive specific explosion energy. We call it resid¬ 
ual recombination energy to make clear that this is energy 
which is not liberated in our simulations, in contrast to the 
energy of th e recombination processes which we identified in 
Section lT^ 



Figure 20. Temporal evolution of the diagnostic energy F'*', the explosion 
energy and the upper limit of the explosion energy also including the 

recombination energy Free for a HC (19.2 M© progenitor) and a LC case 
(20.0 M© progenitor), for PUSH parameters reported in Table[3] 


In Figure |20l we investigate the behavior of the diagnostic 
energy E^{t), and we compare it with our estimate of the ex¬ 


plosion energy £expi(0 = E^^{t) and with its upper limit repre¬ 
sented by E^^ft). We want to emphasize that these quantites 
are obtained from mass integrals above the time-dependent 
mass-c ut, in contrast to most of the energies investigated in 
Section[321 where a fixed mass domain was considered. 

While Fov(0 E^^ft) have already saturated to a con¬ 
stant value at f ^ 1.5 s, even at f « 4.6 s the diagnostic 
energy has not yet converged. £^(0 and E^^ft) approach 
their asymptotic values from below, and any late time increase 
(f > 1.5 s) is due to the energy carried by the neutrino-driven 
wind ejected from the PNS surface. On the other hand, E^{t) 
reaches its maximum around < 1 s after fexpi, when the neu¬ 
trino absorption and the nuclear recombination have released 
most of their energy in the expanding shock wave, and then 
it decreases towards £^(0. since matter with negative total 
specific energy is accreted at the shock. The difference be¬ 
tween £ov(0 £^(0 is mainly represented by the gravi¬ 

tational binding energy of the stellar layers above the shock 
front. Thus, the rate of convergence of the diagnostic energy 
depends on the amount of gravitational binding energy con¬ 
tained in the outer envelope of the star and on the relative 
speed at which the shock propagates inside it. Since the grav¬ 
itational binding energy of the outer layers is similar between 
the two explored models, the different rate of convergence de¬ 
pends mostly on the different expansion velocity of the shock 
w ave, which is larger for more energetic HC models. 

I Yamamoto et alJ (120131) found for a 15 M© progenitor that 
the diagnostic energy saturates and thus reaches the asymp¬ 
totic explosion energy already between 1 and 2 s post-bounce. 
This difference to what we find is related to the different pro¬ 
genitors used and, in particular, to the different binding energy 
of the outer envelopes, which is expected to be much smaller 
for a 15 Mq progenitor t han for a ~ 20 Mp progenitor (see, 
for example. Figure 5 of iBurro'^ l20 1 3l) . Nevertheless, we 
conclude that the diagnostic energy is in general (i.e., without 
further considerations) not suited to give an accurate estimate 
of the explosion energy at early times. 

4.8. Comparison with other works 

A similar fitting to SN 1987A energetics has been done 
for multi-dimensional simulation s (2D and 30) using a light- 
bulb scheme for the neutrinos by iHandv et alJ (l2014h As ini¬ 
tial conditions they used a post-collapse m odel based on th e 
15 Mq blue supergiant progenitor model of lWoosl^ ([198^. 
Even if they did not provide the corresponding compactness, 
the values of the accretion rate (~ 0.2 - 0.3 Mq s“') and of 
the electron neutrino luminosity (~ 1.8 - 3.5 x 10^' erg s“*) at 
the onset of the explosion are more compatible with our FC 
models. In their fitting, only the diagnostic explosion energy 
£^ was used at a time of f nh = 1.5 s when i t is ex pected to 
have saturated to £expi (cf. lYamamoto et al.l (l2013h '). But no 
estimates for the nucleosynthesis yields were given. The time 
when the shock reaches 500 km (which corresponds for us to 
fexpi) is significantly lower in their models (90-140 ms after 
bounce), mainly due to the different extension and evolution 
of the shock during the first 100 ms after core bounce. A 
more detailed quantitative comparison (albeit limited by the 
different dimensionality of the two models) requires to use 
a more similar progenitor. However, the advection timescale 
and the mass in the gain region are larger than the correspond¬ 
ing values we have obtained in all our models, as expected 
from the larger average residence time resulting from multi¬ 
dimensional hydrodynamical effects. 
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lUgliano et alJ (1201 2h also calibrated their spherically sym¬ 
metric exploding models with the observational constraints 
from SN 1987A, and used progenitor mode ls identical to the 
ones we have adopted (IWooslev et al.ll2()()3) . They also found 
that the remnant mass and the properties of the explosion ex¬ 
hibit a large variability inside the narrow 18-21 M© ZAMS 
mass window (they even found some non-exploding models). 
However, they did not find any clear trend with progenitor 
compactness (for example, their calibration model is repre¬ 
sented by the 19.8 M© ZAMS mass progenitor which belongs 
to the LC sample). The explosion timescales for models in 
the 18-21 Mq ZAMS mass interval are much longer in their 
case (fexpi ~ 0.3 - 1 s), while their range for the explosion 
energy (0.6 - 1.6 Bethe) is relatively compatible with ours 
(0.4 - 1.6 Bethe). Clearly, all these differences are related to 
the numerous diversities between the two models. 

A possible relation between explosion properites and 
progenitor compactnes s has been first pointed out by 
lO’ Connor & Ottl (120101) . who searched for a minimum en¬ 
hanced neutrino energy deposition in spherically symmetric 
models. Similarly to us, they found that more compact pro¬ 
genitors require larger heating efficiency to explode. How¬ 
ever, they do not investigate the explosion energy of their 
models. Moreover, they consider it to be unlikely that a model 
which requires tj > 0.23 (^ 2.5 ^ 0.45) will explode in nature. 
In our analysis, we have interpreted a large neutrino heating 
efficiency in spherically symmetric models as an effective way 
to take into account longer residence time inside the gain re¬ 
gion. We have pointed out that HC models, characterized by 
larger //tot, are required to obtain the observed properties of 
SN 1987A. However, these models still have ^ 2.5 < 0.45 and 
our average heat i ng effi ciency are below the critical value of 
lO’Connor&Oa(l2m0h . 

A clear correlation between explosion properties and 
progenitor compactne ss has been recently discussed by 
iNakamura et al.l (12014ft . They performed systematic 2D cal¬ 
culations of exploding CCSNe for a large variety of progen¬ 
itors, using the IDSA to model Ve and Vg transport. Due 
to computational limitations and due to the usage of only a 
NSE EOS, their simulations were limited to ~ 1 s after core 
bounce. Thus, they could not ensure the convergence of the 
diagnostic energy and could not directly compare their results 
with CCSN observables. However, they found trends with 
compactness similar to the ones we have found in our reduce 
sample. 

Other authors have also compared the predicted explosion 
energy and Ni yield from their models to the ob servational 
constraints. Eor example. lYamamoto et al.l (12013ft . using the 
neutrino light-bulb method to trigger explosions in spherical 
symmetry, found a similar trend between explosion energies 
and nickel masses as we found (see Table |6]l. They also com¬ 
pared to a thermal bomb model with similar explosion ener¬ 
gies and mass cut, and found that the neutrino heating mech¬ 
anism leads to systematically larger ^®Ni yields. They re¬ 
lated it to higher peak temperatures, which appear because 
a greater thermal energy is required to unbind the accreting 
envelope. They also concluded that the neutrino-driven mech- 
anis m is more similar to pi ston-driven models by comparing 
with lYoung & Erverl (12007ft . The problem of overproducing 
^^Ni i s lessened in the 2D simulations of lYamamoto et aP 
(12013ft because of slightly lower peak temperatures and the 
occurrence of fallback. 

The conclusions drawn in Section 13.21 about the contribu¬ 


tions of nuclear reactions to the explosion energy are some¬ 
what opposite to what can be found in other w orks in the lit¬ 
erature. Eor example, lYamamoto et al.l (12013ft state that the 
contribution of the nuclear reactions to the explosion energy 
is comparable to or greater than that of neutrino heating. Eur- 
thermore, they identify the recombinations of nucleons into 
nuclei in NSE as the most important nuclear reactions. How¬ 
ever, they also point out that this “recombination energy even¬ 
tually originates from neutrino heating”. We think that this 
aspect is crucial for understanding the global energetics. In¬ 
deed, if we had started the analysis presented in Eigure|8]not 
at bounce but at fexpi we would also have identified a strong 
contribution from the nuclear reactions, given roughly by the 
difference between -(Emass - Emass.o) at fexpi (which is close to 
the minimum) and the final value. However, as is clear from 
the Eigure, roughly the same amount of energy was actually 
taken from the thermal energy before fexpi- The dominant net 
contribution to the explosion energy originates from neutrino 
heating, as is evide nt fro m Eigure|7]and as we haved discussed 
in detail in Section l372l 

5. SUMMARY AND CONCLUSIONS 

The investigation of the explosion mechanism of CCSNe as 
well as accurate explorations of all the aspects related with it, 
is a long lasting, but still fascinating problem. Sophisticated 
multi-dimensional hydrodynamical simulations, possibly in¬ 
cluding detailed neutrino transport, microphysical EOS, mag¬ 
netic fields and aspherical properites of the progenitor struc¬ 
ture, are ultimately required to address this problem. The high 
computational costs of such models and the uncertanties in 
several necessary ingredients still motivate the usage of effec¬ 
tive spherically symmetric models to perform extended pro¬ 
genitor studies. 

In this work we have presented a new method, PUSH, 
for artificially triggering parametrized core-collapse super¬ 
nova explosions of massive stars in spherical symmetry. The 
method provides a robust and computationally affordable 
framework to study important aspects of core-collapse super¬ 
novae that require modeling of the explosion for several sec¬ 
onds after its onset for extended sets of progenitors. Eor ex¬ 
ample, the effects of the shock passage through the star, the 
neutron star mass distribution, the determination of the explo¬ 
sion energy, or explosive supernova nucleosynthesis. Here, 
we have focused on the exploration of basic explosion prop¬ 
erties and on the calibration of PUSH by reproducing observ¬ 
ables of SN 1987A. We considered progenitors in the ZAMS 
mass range of 18 - 21 M 0 which corre sponds to typical values 
for th e progenitor mass of SN 1987A (iShigevama & Nomotol 

dH). 

Unlike traditional methods (such as thermal bombs, pistons, 
or neutrino light-bulbs), our method does not require any ex¬ 
ternal source of energy to trigger the explosion nor a modifi¬ 
cation of the charged-current neutrino reactions. Instead, the 
PUSH method taps a fraction of the energy from muon- and 
tau-neutrinos which are emitted by the PNS. This additional 
energy is deposited inside the gain region for a limited time af¬ 
ter core bounce. The introduction of a local heating term that 
is only active where electron-neutrinos are heating and where 
neutrino-driven convection can occur is inspired by qualita¬ 
tive properties of multi-dimensional CCSN simulations. We 
have two major free parameters, frise, describing the temporal 
evolution of PUSH, and kpush, controlling the strength. They 
are determined by comparing the outcome of our simulations 
with observations. 
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Our setup allows us to model the entire relevant domain, 
including the PNS and the ejecta. In particular, (i) the ther¬ 
modynamic properties of matter both in NSE and non-NSE 
conditions are treated accurately; (ii) the neutrino luminosi¬ 
ties are directly related to the PNS evolution and to the mass 
accretion history; (iii) the evolution of the electron fraction is 
followed by a radiative transport scheme for electron flavor 
neutrinos, which is important for the nucleosynthesis calcula¬ 
tions. 

We have studied the evolution of the explosion energy and 
how it is generated. The energy deposition by neutrinos is the 
main cause of the increase of the total energy of the ejecta 
and, thus, the main source of the explosion energy. The net 
nuclear binding energy released by the ejecta during the whole 
supernova (including both the initial endothermic photodisso¬ 
ciation and the final exothermic explosive burning) is positive, 
but much smaller than the energy provided by neutrinos. Eur- 
thermore, we obtain an approximate convergence of the ex¬ 
plosion energy typically only after 1 to 2 seconds and only if 
the full progenitor structure is taken into account. Vice-versa, 
we And that the so-called “diagnostic energy” is, in general, 
not suited to give an accurate estimate of the explosion energy 
at early times. 

Our broad parameter exploration has revealed a distinction 
between high compactness (^ 1.75 > 0.45) and low compact¬ 
ness (^ 1.75 < 0.45) progenitor models for the ZAMS mass 
range of 18 - 21 Mq. The EC models tend to explode earlier, 
with lower explosion energy, and with a lower remnant mass. 
When the HC models explode, they tend to explode later, 
more energetically, and producing more massive remnants. 
This is due to different accretion histories of the EC and HC 
models. The HC models have larger accretion rates, which 
produce larger neutrino luminosities, (marginally) harder neu¬ 
trino spectra, and a stronger ram pressure at the shock. In or¬ 
der to overcome this pressure a more intense neutrino energy 
deposition is required behind the shock. And, once the ex¬ 
plosion has been launched, a more intense energy deposition 
inside the expanding shock is observed. Thus, HC models re¬ 
quire more time to explode but the resulting explosions are 
more energetic. 

The fitting of the PUSH parameters to observations of 
SN 1987A has lead to several interesting conclusions. The 
requirement of an explosion energy around 1 Bethe has re¬ 
stricted our progenitor search to HC models. At the same 
time, our parameter space exploration has shown that a con¬ 
straint on the explosion energy is equivalent to a tight correla¬ 
tion between the two most relevant PUSH parameters, tnse and 
kpush: if a certain explosion energy has to be achived, a longer 
timescale for PUSH to reach its maximum efficiency (frise) has 
to be compensated by a larger PUSH strength (kpush)- This de¬ 
generacy can be broken by including nucleosynthesis yields in 
the calibration of the free parameters. 

We And an overproduction of ^®Ni for runs with an explo¬ 
sion energy around and above 1 Bethe. This problem is ob¬ 
served for all the tested parameter choices and progenitors that 
provide a sufficiently high explosion energy. Thus, fallback is 
necessary in our models to reproduce the observed nucleosyn¬ 
thesis yields. This fallback could be associated with the for¬ 
mation of a reverse shock when the forward shock reaches the 
hydrogen shell. The relatively large amount of fallback that 
we use (0.1 Mq) is consistent with observational constraints 
from SN 1987A and with explicit calculations of the fallback 
for exploding models of ~ 20 M© ZAMS mass progenitors 


(IChevalieilll989l : IUgliano et alJl2012h . 

The production of sensitive to the electron frac¬ 

tion of the innermost ejecta. A final mass cut initially located 
inside the silicon shell can provide slightly neutron rich ejecta, 
corresponding to the conditions required to fit the 57-58jqj 
yields of SN 1987A. We found that this is only possible for the 
18.0 M© and 19.4 M© ZAMS mass progenitors, whereas for 
the other HC models, characterized by larger ^ 1 . 75 , the mass 
cut is located inside the oxygen shell. The 18.0 M© and 19.4 
M© ZAMS mass progenitors can explain the energetics and all 
nickel yields if fallback is included. Eor ^"^Ti, in contrast, we 
And that it is underproduced. However, we have shown that 
uncertainties in the relevant nuclear reaction rates, together 
with mixing of the ejecta, can help reducing this discrepancy. 

Our work has demonstrated that the progenitor structure 
and composition are of great importance for the nucleosynthe¬ 
sis yields. Recently, it has been pointed out that asphericities 
in the progenitor structure could aid th e multi-dimensional 
neutrino-driven supernova mechanism (ICouch & OttI 120131 : 
iMuller & Jankall2015l ICouch et alJl2015l) . Eor our work, the 
compositional changes indu ced by multi-dimen sional effects 
in the progenitor evolution (lArnett et al.ll2015h would be of 
particular interest and could be the subject of future work. 
However, at present, databases with large sets of progenitors 
are only available for calculations that were done in spherical 
symmetry. 

Einally, we have identifled a progenitor (18.0 M© ZAMS 
mass, compactness fi .75 = 0.463 at collapse) that flts the ob¬ 
servables of SN 1987A for a suitable choice of the PUSH pa¬ 
rameters (ton = 80 ms, frise = 200 ms, and kpush = 3.5) and 
assuming 0.1 M© of fallback. The associated explosion en¬ 
ergy is Eexpi = 1.092 Bethe, while M{ ^®Ni) = 0.073 M©. The 
formation of a BH in SN 1987A seems to be unlikely, since 
it would require a much larger fallback compared with our 
analysis and/or an extremely asymmetric explosion. Instead, 
we predict that in SN 1987A a neutron star with a baryonic 
mass of 1 .66 M© was born, corresponding to a gravitational 
mass of 1.50 M© for a cold neutron star with our choice of 
the EOS. This will h opefully be probed by observations soon 
(iZanardo et alJl20T^ . 

Eor our best model of SN 1987A the explosion happens 
on a timescale of a few hundereds of milliseconds after core 
bounce. This timescale is consistent with the overall picture 
of a neutrino-driven supernova, and broadly compatible with 
the Arst results obtained in exploding, self-consistent, multi¬ 
dimensional simulations. 

Erom exploring the progenitor range of 18 - 21 M© ZAMS 
mass we found indications of a correlation between explo¬ 
sion properties and the compactness of the progenitor model. 
However, a more complete analysis will require the explo¬ 
ration of a larger set of progenitors with the PUSH method. 
This will be the topic of a forthcoming work. An extended 
study considering all possible progenitors for core-collapse 
supemovae in the mass range of 8 - 100 M© will be rele¬ 
vant for several open questions in nuclear astrophysics, for 
example for the comparison of predicted to observed explo¬ 
sion ener gies, neutron-star rem nant masses, and ejected ^®Ni 
(see, e.g.. lBruenn et al.l (120141) 1 or for the prediction of com¬ 
plete nucleosynthesis yields of all elements which is a crucial 
input to galactic chemical evolution. A full progenitor study 
could also give more insight into the extent to which the com¬ 
pactness parameter affects the supernova energetics and nu¬ 
cleosynthesis. 
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